From eeae3ec78365fd7a1f1ad661c100e1108401b0cd Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 3 Mar 2014 12:23:20 -0800 Subject: [PATCH] Changed LogicallyOrthogonalMesh to LogicallyRectMesh and updated all dependencies. LOM --> LRM removed LomView.py, and put plot grid code inside Mesh code. Added tutorial style introduction to the mesh. --- SimPEG/Mesh/InnerProducts.py | 21 +- ...OrthogonalMesh.py => LogicallyRectMesh.py} | 119 ++++++++++- SimPEG/Mesh/LomView.py | 104 ---------- SimPEG/Mesh/__init__.py | 3 +- SimPEG/Tests/TestUtils.py | 12 +- SimPEG/Tests/test_LogicallyOrthogonalMesh.py | 104 ---------- SimPEG/Tests/test_LogicallyRectMesh.py | 104 ++++++++++ SimPEG/Tests/test_massMatrices.py | 4 +- SimPEG/Tests/test_operators.py | 6 +- SimPEG/Utils/__init__.py | 4 +- SimPEG/Utils/{lomutils.py => lrmutils.py} | 0 SimPEG/Utils/meshutils.py | 2 +- docs/api_Mesh.rst | 192 +++++++++++++++--- docs/api_MeshCode.rst | 35 ++++ docs/api_Utils.rst | 4 +- 15 files changed, 445 insertions(+), 269 deletions(-) rename SimPEG/Mesh/{LogicallyOrthogonalMesh.py => LogicallyRectMesh.py} (74%) delete mode 100644 SimPEG/Mesh/LomView.py delete mode 100644 SimPEG/Tests/test_LogicallyOrthogonalMesh.py create mode 100644 SimPEG/Tests/test_LogicallyRectMesh.py rename SimPEG/Utils/{lomutils.py => lrmutils.py} (100%) create mode 100644 docs/api_MeshCode.rst diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 1c85e406..99490674 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -16,6 +16,7 @@ class InnerProducts(object): :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :param bool returnP: returns the projection matrices :param bool invertProperty: inverts the material property + :param bool doFast: do a faster implementation if available. :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ @@ -93,11 +94,13 @@ class InnerProducts(object): return self._getInnerProductDeriv(materialProperty, v, P, self.nF) - def getEdgeInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False, doFast=True): + def getEdgeInnerProduct(self, materialProperty=None, returnP=False, + invertProperty=False, doFast=True): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :param bool returnP: returns the projection matrices :param bool invertProperty: inverts the material property + :param bool doFast: do a faster implementation if available. :rtype: scipy.csr_matrix :return: M, the inner product matrix (nE, nE) """ @@ -353,7 +356,7 @@ def _getFacePxx_Rectangular(M): iijj = ndgrid(i, j) ii, jj = iijj[:, 0], iijj[:, 1] - if M._meshType == 'LOM': + if M._meshType == 'LRM': fN1 = M.r(M.normals, 'F', 'Fx', 'M') fN2 = M.r(M.normals, 'F', 'Fy', 'M') @@ -378,7 +381,7 @@ def _getFacePxx_Rectangular(M): PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF)) - if M._meshType == 'LOM': + if M._meshType == 'LRM': I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]), getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy])) PXX = I2x2 * PXX @@ -401,7 +404,7 @@ def _getFacePxxx_Rectangular(M): iijjkk = ndgrid(i, j, k) ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - if M._meshType == 'LOM': + if M._meshType == 'LRM': fN1 = M.r(M.normals, 'F', 'Fx', 'M') fN2 = M.r(M.normals, 'F', 'Fy', 'M') fN3 = M.r(M.normals, 'F', 'Fz', 'M') @@ -435,7 +438,7 @@ def _getFacePxxx_Rectangular(M): PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr() - if M._meshType == 'LOM': + if M._meshType == 'LRM': I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]), getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]), getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ])) @@ -450,7 +453,7 @@ def _getEdgePxx_Rectangular(M): iijj = ndgrid(i, j) ii, jj = iijj[:, 0], iijj[:, 1] - if M._meshType == 'LOM': + if M._meshType == 'LRM': eT1 = M.r(M.tangents, 'E', 'Ex', 'M') eT2 = M.r(M.tangents, 'E', 'Ey', 'M') @@ -470,7 +473,7 @@ def _getEdgePxx_Rectangular(M): PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr() - if M._meshType == 'LOM': + if M._meshType == 'LRM': I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]), getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j])) PXX = I2x2 * PXX @@ -484,7 +487,7 @@ def _getEdgePxxx_Rectangular(M): iijjkk = ndgrid(i, j, k) ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - if M._meshType == 'LOM': + if M._meshType == 'LRM': eT1 = M.r(M.tangents, 'E', 'Ex', 'M') eT2 = M.r(M.tangents, 'E', 'Ey', 'M') eT3 = M.r(M.tangents, 'E', 'Ez', 'M') @@ -513,7 +516,7 @@ def _getEdgePxxx_Rectangular(M): PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr() - if M._meshType == 'LOM': + if M._meshType == 'LRM': I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]), getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]), getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k])) diff --git a/SimPEG/Mesh/LogicallyOrthogonalMesh.py b/SimPEG/Mesh/LogicallyRectMesh.py similarity index 74% rename from SimPEG/Mesh/LogicallyOrthogonalMesh.py rename to SimPEG/Mesh/LogicallyRectMesh.py index 7b4ccd73..d606fe78 100644 --- a/SimPEG/Mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/Mesh/LogicallyRectMesh.py @@ -2,7 +2,6 @@ from SimPEG import Utils, np from BaseMesh import BaseRectangularMesh from DiffOperators import DiffOperators from InnerProducts import InnerProducts -from LomView import LomView # Some helper functions. length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5 @@ -11,24 +10,24 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2)) normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2)) -class LogicallyOrthogonalMesh(BaseRectangularMesh, DiffOperators, InnerProducts, LomView): +class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts): """ - LogicallyOrthogonalMesh is a mesh class that deals with logically orthogonal meshes. + LogicallyRectMesh is a mesh class that deals with logically rectangular meshes. - Example of a logically orthogonal mesh: + Example of a logically rectangular mesh: .. plot:: :include-source: from SimPEG import Mesh, Utils - X, Y = Utils.exampleLomGird([3,3],'rotate') - M = Mesh.LogicallyOrthogonalMesh([X, Y]) + X, Y = Utils.exampleLrmGrid([3,3],'rotate') + M = Mesh.LogicallyRectMesh([X, Y]) M.plotGrid(showIt=True) """ __metaclass__ = Utils.SimPEGMetaClass - _meshType = 'LOM' + _meshType = 'LRM' def __init__(self, nodes): assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray" @@ -39,7 +38,7 @@ class LogicallyOrthogonalMesh(BaseRectangularMesh, DiffOperators, InnerProducts, assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % i) assert len(nodes[0].shape) == len(nodes), "Dimension mismatch" - assert len(nodes[0].shape) > 1, "Not worth using LOM for a 1D mesh." + assert len(nodes[0].shape) > 1, "Not worth using LRM for a 1D mesh." BaseRectangularMesh.__init__(self, np.array(nodes[0].shape)-1, None) @@ -329,6 +328,106 @@ class LogicallyOrthogonalMesh(BaseRectangularMesh, DiffOperators, InnerProducts, _tangents = None tangents = property(**tangents()) + + + ############################################# + # Plotting Functions # + ############################################# + + def plotGrid(self, length=0.05, showIt=False): + """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. + + + .. plot:: + :include-source: + + from SimPEG import Mesh, Utils + X, Y = Utils.exampleLrmGrid([3,3],'rotate') + M = Mesh.LogicallyRectMesh([X, Y]) + M.plotGrid(showIt=True) + + """ + import matplotlib.pyplot as plt + import matplotlib + from mpl_toolkits.mplot3d import Axes3D + mkvc = Utils.mkvc + + 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() + + X = np.r_[X1, X2] + Y = np.r_[Y1, Y2] + + plt.plot(X, Y) + + 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') + + plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo') + + 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-') + + 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-') + + 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.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') + + 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() + + 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() + Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten() + + X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten() + Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten() + Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*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.set_zlabel('x3') + + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + + if showIt: plt.show() + + if __name__ == '__main__': nc = 5 h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)]) @@ -338,9 +437,9 @@ if __name__ == '__main__': dee3 = True if dee3: X, Y, Z = Utils.ndgrid(h1, h2, h3, vector=False) - M = LogicallyOrthogonalMesh([X, Y, Z]) + M = LogicallyRectMesh([X, Y, Z]) else: X, Y = Utils.ndgrid(h1, h2, vector=False) - M = LogicallyOrthogonalMesh([X, Y]) + M = LogicallyRectMesh([X, Y]) print M.r(M.normals, 'F', 'Fx', 'V') diff --git a/SimPEG/Mesh/LomView.py b/SimPEG/Mesh/LomView.py deleted file mode 100644 index 4eceb918..00000000 --- a/SimPEG/Mesh/LomView.py +++ /dev/null @@ -1,104 +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 LomView(object): - """ - Provides viewing functions for LogicallyOrthogonalMesh - - This class is inherited by LogicallyOrthogonalMesh - - """ - def __init__(self): - pass - - def plotGrid(self, length=0.05, showIt=False): - """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. - - - .. plot:: - :include-source: - - from SimPEG import Mesh, Utils - X, Y = Utils.exampleLomGird([3,3],'rotate') - M = Mesh.LogicallyOrthogonalMesh([X, Y]) - M.plotGrid(showIt=True) - - """ - 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() - - X = np.r_[X1, X2] - Y = np.r_[Y1, Y2] - - plt.plot(X, Y) - - 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') - - plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo') - - 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-') - - 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-') - - 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.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') - - 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() - - 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() - Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten() - - X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten() - Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten() - Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*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.set_zlabel('x3') - - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - ax.set_ylabel('x2') - - if showIt: plt.show() diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py index 47b30243..0fa60201 100644 --- a/SimPEG/Mesh/__init__.py +++ b/SimPEG/Mesh/__init__.py @@ -1,9 +1,8 @@ from Cyl1DMesh import Cyl1DMesh from TensorMesh import TensorMesh from TreeMesh import TreeMesh -from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh +from LogicallyRectMesh import LogicallyRectMesh from BaseMesh import BaseMesh, BaseRectangularMesh from TensorView import TensorView -from LomView import LomView from InnerProducts import InnerProducts from DiffOperators import DiffOperators diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index b46fbf10..da18f537 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt from numpy.linalg import norm from SimPEG.Utils import mkvc, sdiag from SimPEG import Utils -from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh +from SimPEG.Mesh import TensorMesh, LogicallyRectMesh import numpy as np import scipy.sparse as sp import unittest @@ -104,7 +104,7 @@ class OrderTest(unittest.TestCase): max_h = max([np.max(hi) for hi in self.M.h]) return max_h - elif 'LOM' in self._meshType: + elif 'LRM' in self._meshType: if 'uniform' in self._meshType: kwrd = 'rect' elif 'rotate' in self._meshType: @@ -114,11 +114,11 @@ class OrderTest(unittest.TestCase): if self.meshDimension == 1: raise Exception('Lom not supported for 1D') elif self.meshDimension == 2: - X, Y = Utils.exampleLomGird([nc, nc], kwrd) - self.M = LogicallyOrthogonalMesh([X, Y]) + X, Y = Utils.exampleLrmGrid([nc, nc], kwrd) + self.M = LogicallyRectMesh([X, Y]) elif self.meshDimension == 3: - X, Y, Z = Utils.exampleLomGird([nc, nc, nc], kwrd) - self.M = LogicallyOrthogonalMesh([X, Y, Z]) + X, Y, Z = Utils.exampleLrmGrid([nc, nc, nc], kwrd) + self.M = LogicallyRectMesh([X, Y, Z]) return 1./nc def getError(self): diff --git a/SimPEG/Tests/test_LogicallyOrthogonalMesh.py b/SimPEG/Tests/test_LogicallyOrthogonalMesh.py deleted file mode 100644 index 0c5b3247..00000000 --- a/SimPEG/Tests/test_LogicallyOrthogonalMesh.py +++ /dev/null @@ -1,104 +0,0 @@ -import numpy as np -import unittest -from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh -from SimPEG.Utils import ndgrid - - -class BasicLOMTests(unittest.TestCase): - - def setUp(self): - a = np.array([1, 1, 1]) - b = np.array([1, 2]) - c = np.array([1, 4]) - gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h] - X, Y = ndgrid(gridIt([a, b]), vector=False) - self.TM2 = TensorMesh([a, b]) - self.LOM2 = LogicallyOrthogonalMesh([X, Y]) - X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False) - self.TM3 = TensorMesh([a, b, c]) - self.LOM3 = LogicallyOrthogonalMesh([X, Y, Z]) - - def test_area_3D(self): - test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2]) - self.assertTrue(np.all(self.LOM3.area == test_area)) - - def test_vol_3D(self): - test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8]) - np.testing.assert_almost_equal(self.LOM3.vol, test_vol) - self.assertTrue(True) # Pass if you get past the assertion. - - def test_vol_2D(self): - test_vol = np.array([1, 1, 1, 2, 2, 2]) - t1 = np.all(self.LOM2.vol == test_vol) - self.assertTrue(t1) - - def test_edge_3D(self): - test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]) - t1 = np.all(self.LOM3.edge == test_edge) - self.assertTrue(t1) - - def test_edge_2D(self): - test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2]) - t1 = np.all(self.LOM2.edge == test_edge) - self.assertTrue(t1) - - def test_tangents(self): - T = self.LOM2.tangents - self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM2.nEx))) - self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM2.nEx))) - self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM2.nEy))) - self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM2.nEy))) - - T = self.LOM3.tangents - self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM3.nEx))) - self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM3.nEx))) - self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LOM3.nEx))) - - self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM3.nEy))) - self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM3.nEy))) - self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LOM3.nEy))) - - self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LOM3.nEz))) - self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LOM3.nEz))) - self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LOM3.nEz))) - - def test_normals(self): - N = self.LOM2.normals - self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM2.nFx))) - self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM2.nFx))) - self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM2.nFy))) - self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM2.nFy))) - - N = self.LOM3.normals - self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM3.nFx))) - self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM3.nFx))) - self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LOM3.nFx))) - - self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM3.nFy))) - self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM3.nFy))) - self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LOM3.nFy))) - - self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LOM3.nFz))) - self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LOM3.nFz))) - self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LOM3.nFz))) - - def test_grid(self): - self.assertTrue(np.all(self.LOM2.gridCC == self.TM2.gridCC)) - self.assertTrue(np.all(self.LOM2.gridN == self.TM2.gridN)) - self.assertTrue(np.all(self.LOM2.gridFx == self.TM2.gridFx)) - self.assertTrue(np.all(self.LOM2.gridFy == self.TM2.gridFy)) - self.assertTrue(np.all(self.LOM2.gridEx == self.TM2.gridEx)) - self.assertTrue(np.all(self.LOM2.gridEy == self.TM2.gridEy)) - - self.assertTrue(np.all(self.LOM3.gridCC == self.TM3.gridCC)) - self.assertTrue(np.all(self.LOM3.gridN == self.TM3.gridN)) - self.assertTrue(np.all(self.LOM3.gridFx == self.TM3.gridFx)) - self.assertTrue(np.all(self.LOM3.gridFy == self.TM3.gridFy)) - self.assertTrue(np.all(self.LOM3.gridFz == self.TM3.gridFz)) - self.assertTrue(np.all(self.LOM3.gridEx == self.TM3.gridEx)) - self.assertTrue(np.all(self.LOM3.gridEy == self.TM3.gridEy)) - self.assertTrue(np.all(self.LOM3.gridEz == self.TM3.gridEz)) - - -if __name__ == '__main__': - unittest.main() diff --git a/SimPEG/Tests/test_LogicallyRectMesh.py b/SimPEG/Tests/test_LogicallyRectMesh.py new file mode 100644 index 00000000..5d223f68 --- /dev/null +++ b/SimPEG/Tests/test_LogicallyRectMesh.py @@ -0,0 +1,104 @@ +import numpy as np +import unittest +from SimPEG.Mesh import TensorMesh, LogicallyRectMesh +from SimPEG.Utils import ndgrid + + +class BasicLRMTests(unittest.TestCase): + + def setUp(self): + a = np.array([1, 1, 1]) + b = np.array([1, 2]) + c = np.array([1, 4]) + gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h] + X, Y = ndgrid(gridIt([a, b]), vector=False) + self.TM2 = TensorMesh([a, b]) + self.LRM2 = LogicallyRectMesh([X, Y]) + X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False) + self.TM3 = TensorMesh([a, b, c]) + self.LRM3 = LogicallyRectMesh([X, Y, Z]) + + def test_area_3D(self): + test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2]) + self.assertTrue(np.all(self.LRM3.area == test_area)) + + def test_vol_3D(self): + test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8]) + np.testing.assert_almost_equal(self.LRM3.vol, test_vol) + self.assertTrue(True) # Pass if you get past the assertion. + + def test_vol_2D(self): + test_vol = np.array([1, 1, 1, 2, 2, 2]) + t1 = np.all(self.LRM2.vol == test_vol) + self.assertTrue(t1) + + def test_edge_3D(self): + test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]) + t1 = np.all(self.LRM3.edge == test_edge) + self.assertTrue(t1) + + def test_edge_2D(self): + test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2]) + t1 = np.all(self.LRM2.edge == test_edge) + self.assertTrue(t1) + + def test_tangents(self): + T = self.LRM2.tangents + self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM2.nEx))) + self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM2.nEx))) + self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM2.nEy))) + self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM2.nEy))) + + T = self.LRM3.tangents + self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM3.nEx))) + self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM3.nEx))) + self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LRM3.nEx))) + + self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM3.nEy))) + self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM3.nEy))) + self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LRM3.nEy))) + + self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LRM3.nEz))) + self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LRM3.nEz))) + self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LRM3.nEz))) + + def test_normals(self): + N = self.LRM2.normals + self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM2.nFx))) + self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM2.nFx))) + self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM2.nFy))) + self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM2.nFy))) + + N = self.LRM3.normals + self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM3.nFx))) + self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM3.nFx))) + self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LRM3.nFx))) + + self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM3.nFy))) + self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM3.nFy))) + self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LRM3.nFy))) + + self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LRM3.nFz))) + self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LRM3.nFz))) + self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LRM3.nFz))) + + def test_grid(self): + self.assertTrue(np.all(self.LRM2.gridCC == self.TM2.gridCC)) + self.assertTrue(np.all(self.LRM2.gridN == self.TM2.gridN)) + self.assertTrue(np.all(self.LRM2.gridFx == self.TM2.gridFx)) + self.assertTrue(np.all(self.LRM2.gridFy == self.TM2.gridFy)) + self.assertTrue(np.all(self.LRM2.gridEx == self.TM2.gridEx)) + self.assertTrue(np.all(self.LRM2.gridEy == self.TM2.gridEy)) + + self.assertTrue(np.all(self.LRM3.gridCC == self.TM3.gridCC)) + self.assertTrue(np.all(self.LRM3.gridN == self.TM3.gridN)) + self.assertTrue(np.all(self.LRM3.gridFx == self.TM3.gridFx)) + self.assertTrue(np.all(self.LRM3.gridFy == self.TM3.gridFy)) + self.assertTrue(np.all(self.LRM3.gridFz == self.TM3.gridFz)) + self.assertTrue(np.all(self.LRM3.gridEx == self.TM3.gridEx)) + self.assertTrue(np.all(self.LRM3.gridEy == self.TM3.gridEy)) + self.assertTrue(np.all(self.LRM3.gridEz == self.TM3.gridEz)) + + +if __name__ == '__main__': + unittest.main() diff --git a/SimPEG/Tests/test_massMatrices.py b/SimPEG/Tests/test_massMatrices.py index 200d2884..3d3946ad 100644 --- a/SimPEG/Tests/test_massMatrices.py +++ b/SimPEG/Tests/test_massMatrices.py @@ -6,7 +6,7 @@ from TestUtils import OrderTest class TestInnerProducts(OrderTest): """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" - meshTypes = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM'] + meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM'] meshDimension = 3 meshSizes = [16, 32] @@ -97,7 +97,7 @@ class TestInnerProducts(OrderTest): class TestInnerProducts2D(OrderTest): """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" - meshTypes = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM'] + meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM'] meshDimension = 2 meshSizes = [4, 8, 16, 32, 64, 128] diff --git a/SimPEG/Tests/test_operators.py b/SimPEG/Tests/test_operators.py index 67741aa7..1cf99402 100644 --- a/SimPEG/Tests/test_operators.py +++ b/SimPEG/Tests/test_operators.py @@ -3,7 +3,7 @@ import unittest from TestUtils import OrderTest import matplotlib.pyplot as plt -MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM'] +MESHTYPES = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM'] call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] @@ -37,7 +37,7 @@ class TestCurl(OrderTest): curlE_anal = self.M.projectFaceVector(Fc) curlE = self.M.edgeCurl.dot(E) - if self._meshType == 'rotateLOM': + if self._meshType == 'rotateLRM': # Really it is the integration we should be caring about: # So, let us look at the l2 norm. err = np.linalg.norm(self.M.area*(curlE - curlE_anal), 2) @@ -207,7 +207,7 @@ class TestFaceDiv3D(OrderTest): divF = self.M.faceDiv.dot(F) divF_anal = call3(sol, self.M.gridCC) - if self._meshType == 'rotateLOM': + if self._meshType == 'rotateLRM': # Really it is the integration we should be caring about: # So, let us look at the l2 norm. err = np.linalg.norm(self.M.vol*(divF-divF_anal), 2) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index f45a1039..a44cfd74 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -1,6 +1,6 @@ from matutils import * -from meshutils import exampleLomGird, meshTensors -from lomutils import volTetra, faceInfo, indexCube +from meshutils import exampleLrmGrid, meshTensors +from lrmutils import volTetra, faceInfo, indexCube from interputils import interpmat from ipythonutils import easyAnimate as animate import ModelBuilder diff --git a/SimPEG/Utils/lomutils.py b/SimPEG/Utils/lrmutils.py similarity index 100% rename from SimPEG/Utils/lomutils.py rename to SimPEG/Utils/lrmutils.py diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index f8add82b..48dcb8cf 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -2,7 +2,7 @@ import numpy as np from scipy import sparse as sp from matutils import mkvc, ndgrid, sub2ind, sdiag -def exampleLomGird(nC, exType): +def exampleLrmGrid(nC, exType): assert type(nC) == list, "nC must be a list containing the number of nodes" assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions" exType = exType.lower() diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst index 2a3443d7..46bd877b 100644 --- a/docs/api_Mesh.rst +++ b/docs/api_Mesh.rst @@ -1,38 +1,182 @@ .. _api_Mesh: +SimPEG Meshes +************* -Tensor Mesh -=========== +The Mesh objects in SimPEG provide a numerical grid on which to solve +differential equations. Each mesh type has a similar API to make switching +between different meshes relatively simple. -.. automodule:: SimPEG.Mesh.TensorMesh - :show-inheritance: - :members: - :undoc-members: +Overview of Meshes Available +============================ + +The following meshes are available for use: + +.. toctree:: + :maxdepth: 2 + + api_MeshCode + +Each mesh code follows the guiding principles that are present in this +tutorial, but the details, advantages and disadvantages differ between +the implementations. -Cylindrical 1D Mesh -=================== - -.. automodule:: SimPEG.Mesh.Cyl1DMesh - :show-inheritance: - :members: - :undoc-members: -Logically Orthogonal Mesh -========================= +Variable Locations and Terminology +================================== -.. automodule:: SimPEG.Mesh.LogicallyOrthogonalMesh - :show-inheritance: - :members: - :undoc-members: +We will go over the basics of using a TensorMesh, but these skills are transferable +to the other meshes available in SimPEG. All of the mesh generation code is located +in the Mesh package in SimPEG (i.e. SimPEG.Mesh). -Base Mesh -========= +To create a TensorMesh we need to create mesh tensors, the widths of +each cell of the mesh in each dimension. We will call these tensors h, +and these will be define the constant widths of cells in each dimension +of the TensorMesh. -.. automodule:: SimPEG.Mesh.BaseMesh - :members: - :undoc-members: +.. plot:: + :include-source: + + from SimPEG import Mesh, np + hx = np.r_[3,2,1,1,1,1,2,3] + hy = np.r_[3,1,1,3] + M = Mesh.TensorMesh([hx, hy]) + M.plotGrid(centers=True) +In this simple mesh, the hx vector defines the widths of the cell +in the x dimension, and starts counting from the origin (0,0). The +resulting mesh is divided into cells, and the cell-centers are +plotted above as red circles. Other terminology for this mesh are: + +- cell-centers +- nodes +- faces +- edges + +.. plot:: + :include-source: + + from SimPEG import Mesh, np + import matplotlib.pyplot as plt + hx = np.r_[3,2,1,1,1,1,2,3] + hy = np.r_[3,1,1,3] + M = Mesh.TensorMesh([hx, hy]) + M.plotGrid(faces=True, nodes=True) + plt.title('Cell faces in the x- and y-directions.') + plt.legend(('Nodes', 'X-Faces', 'Y-Faces')) + +Generally, the faces are used to discretize fluxes, quantities that +leave or enter the cells. As such, these fluxes have a direction to +them, which is normal to the cell (i.e. directly out of the cell face). +The plot above shows that x-faces point in the x-direction, and +y-faces point in the y-direction. The nodes are shown in blue, +and lie at the intersection of the grid lines. In a two-dimensional +mesh, the edges actually live in the same location as the faces, +however, they align (or are tangent to) the face. This is easier to +see in 3D, when the edges do not live in the same location as the faces. +In the 3D plot below, the edge variables are seen as black triangles, +and live on the edges(!) of the cell. + +.. plot:: + :include-source: + + from SimPEG import Mesh + Mesh.TensorMesh([1,1,1]).plotGrid(faces=True, edges=True, centers=True) + +How many of each? +----------------- + +When making variables that live in each of these locations, it is +important to know how many of each variable type you are dealing with. +SimPEG makes this pretty easy: + +:: + + In [1]: print M + ---- 2-D TensorMesh ---- + x0: 0.00 + y0: 0.00 + nCx: 8 + nCy: 4 + hx: 3.00, 2.00, 4*1.00, 2.00, 3.00 + hy: 3.00, 2*1.00, 3.00 + + In [2]: count = {'numCells': M.nC, + ....: 'numCells_xDir': M.nCx, + ....: 'numCells_yDir': M.nCy, + ....: 'numCells_vector': M.vnC} + + In [3]: print 'This mesh has %(numCells)d cells, which is %(numCells_xDir)d*%(numCells_yDir)d!!' % count + + This mesh has 32 cells, which is 8*4!! + + In [4]: print count + + { + 'numCells_vector': array([8, 4]), + 'numCells_yDir': 4, + 'numCells_xDir': 8, + 'numCells': 32 + } + +SimPEG also counts the nodes, faces, and edges. + +:: + + Nodes: M.nN, M.nNx, M.nNy, M.nNz, M.vnN + Faces: M.nF, M.nFx, M.nFy, M.nFz, M.vnF, M.vnFx, M.vnFy, M.vnFz + Edges: M.nE, M.nEx, M.nEy, M.nEz, M.vnE, M.vnEx, M.vnEy, M.vnEz + +Face and edge variables have different counts depending on +the dimension of the direction that you are interested in. +In a 4x5 mesh, for example, there is a 5x5 grid of x-faces, +and a 4x6 grid of y-faces. You can count them below! +As such, the vnF(x,y,z) and vnE(x,y,z) properties give the +vector grid size. + +.. plot:: + :include-source: + + from SimPEG import Mesh + Mesh.TensorMesh([4,5]).plotGrid(faces=True) + + +Making Tensors +-------------- + +For tensor meshes, there are some additional functions that can come +in handy. For example, creating mesh tensors can be a bit time +consuming, these can be created speedily by just giving numbers +and sizes of padding. See the example below, that follows this +notation:: + + h1 = ( + (numPad, sizeStart [, increaseFactor]), + (numCore, sizeCode), + (numPad, sizeStart [, increaseFactor]) + ) + +.. plot:: + :include-source: + + from SimPEG import Mesh, Utils + h1 = (5, 10, 1.5), (20, 5), (3, 10) + M = Mesh.TensorMesh(Utils.meshTensors(h1, h1)) + M.plotGrid() + +Hopefully, you now know how to create TensorMesh objects in SimPEG, +and by extension you are also familiar with how to create and use +other types of meshes in this SimPEG framework. + + +The API +======= + +.. toctree:: + :maxdepth: 2 + + api_MeshCode diff --git a/docs/api_MeshCode.rst b/docs/api_MeshCode.rst new file mode 100644 index 00000000..090bb6c2 --- /dev/null +++ b/docs/api_MeshCode.rst @@ -0,0 +1,35 @@ +.. _api_MeshCode: + +Tensor Mesh +=========== + +.. automodule:: SimPEG.Mesh.TensorMesh + :show-inheritance: + :members: + :undoc-members: + + +Cylindrical 1D Mesh +=================== + +.. automodule:: SimPEG.Mesh.Cyl1DMesh + :show-inheritance: + :members: + :undoc-members: + + +Logically Rectangular Mesh +========================== + +.. automodule:: SimPEG.Mesh.LogicallyRectMesh + :show-inheritance: + :members: + :undoc-members: + + +Base Mesh +========= + +.. automodule:: SimPEG.Mesh.BaseMesh + :members: + :undoc-members: diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst index c8dfe47b..f9e24e04 100644 --- a/docs/api_Utils.rst +++ b/docs/api_Utils.rst @@ -15,10 +15,10 @@ Matrix Utilities :members: :undoc-members: -LOM Utilities +LRM Utilities ============= -.. automodule:: SimPEG.Utils.lomutils +.. automodule:: SimPEG.Utils.lrmutils :members: :undoc-members: