From 44a846cb23621a6409d0d4fd612168e5ca7f3b8c Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 30 Sep 2013 11:36:09 -0700 Subject: [PATCH 1/4] added lines to the tensor view. --- SimPEG/TensorView.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/SimPEG/TensorView.py b/SimPEG/TensorView.py index 2b02bfaa..168986a0 100644 --- a/SimPEG/TensorView.py +++ b/SimPEG/TensorView.py @@ -2,6 +2,7 @@ import numpy as np import matplotlib.pyplot as plt import matplotlib from mpl_toolkits.mplot3d import Axes3D +from utils import mkvc class TensorView(object): @@ -213,6 +214,17 @@ class TensorView(object): ax.plot(xc[:, 0], xc[:, 1], 'ro') ax.plot(xs1[:, 0], xs1[:, 1], 'g>') ax.plot(xs2[:, 0], xs2[:, 1], 'g^') + + # Plot the grid 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') @@ -241,6 +253,23 @@ class TensorView(object): 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 + 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') From dde447423493cfe479ef2e06f1f029a1c284bf5a Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 30 Sep 2013 14:53:22 -0700 Subject: [PATCH 2/4] updated plots and documentation --- SimPEG/LogicallyOrthogonalMesh.py | 4 + SimPEG/TensorMesh.py | 5 +- SimPEG/TensorView.py | 82 +++++++++++-------- .../mesh/plot_LogicallyOrthogonalMesh.py | 6 ++ docs/examples/mesh/plot_TensorMesh.py | 21 +++++ docs/examples/mesh/plot_grid_2D.py | 2 +- docs/examples/mesh/plot_grid_3D.py | 2 +- 7 files changed, 83 insertions(+), 39 deletions(-) create mode 100644 docs/examples/mesh/plot_LogicallyOrthogonalMesh.py create mode 100644 docs/examples/mesh/plot_TensorMesh.py diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/LogicallyOrthogonalMesh.py index cfa8b7ed..49e2ee07 100644 --- a/SimPEG/LogicallyOrthogonalMesh.py +++ b/SimPEG/LogicallyOrthogonalMesh.py @@ -16,6 +16,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): """ LogicallyOrthogonalMesh is a mesh class that deals with logically orthogonal meshes. + Example of a logically orthogonal mesh: + + .. plot:: examples/mesh/plot_LogicallyOrthogonalMesh.py + """ _meshType = 'LOM' diff --git a/SimPEG/TensorMesh.py b/SimPEG/TensorMesh.py index c59f4acd..2015838c 100644 --- a/SimPEG/TensorMesh.py +++ b/SimPEG/TensorMesh.py @@ -21,8 +21,9 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): mesh = TensorMesh([hx, hy, hz]) - .. math:: - x^2 = 5 + Example of a padded tensor mesh: + + .. plot:: examples/mesh/plot_TensorMesh.py """ _meshType = 'TENSOR' diff --git a/SimPEG/TensorView.py b/SimPEG/TensorView.py index 168986a0..e0c570cd 100644 --- a/SimPEG/TensorView.py +++ b/SimPEG/TensorView.py @@ -177,9 +177,16 @@ class TensorView(object): if showIt: plt.show() return ph - def plotGrid(self, showIt=False): + 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: @@ -210,20 +217,22 @@ class TensorView(object): xs2 = self.gridFy ax.hold(True) - ax.plot(xn[:, 0], xn[:, 1], 'bs') - ax.plot(xc[:, 0], xc[:, 1], 'ro') - ax.plot(xs1[:, 0], xs1[:, 1], 'g>') - ax.plot(xs2[:, 0], xs2[:, 1], 'g^') + 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 - 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) + 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) @@ -245,30 +254,33 @@ class TensorView(object): xes3 = self.gridEz ax.hold(True) - ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2]) - ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2]) - 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(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]) + 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 - 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) + 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) diff --git a/docs/examples/mesh/plot_LogicallyOrthogonalMesh.py b/docs/examples/mesh/plot_LogicallyOrthogonalMesh.py new file mode 100644 index 00000000..55350946 --- /dev/null +++ b/docs/examples/mesh/plot_LogicallyOrthogonalMesh.py @@ -0,0 +1,6 @@ +from SimPEG import LogicallyOrthogonalMesh, utils +import matplotlib.pyplot as plt +X, Y = utils.exampleLomGird([3,3],'rotate') +M = LogicallyOrthogonalMesh([X, Y]) +M.plotGrid() +plt.show() diff --git a/docs/examples/mesh/plot_TensorMesh.py b/docs/examples/mesh/plot_TensorMesh.py new file mode 100644 index 00000000..94c627f8 --- /dev/null +++ b/docs/examples/mesh/plot_TensorMesh.py @@ -0,0 +1,21 @@ +import numpy as np +import matplotlib.pyplot as plt +from SimPEG import TensorMesh + +pad = 7 +padfactor = 1.4 +xpad = (np.ones(pad)*padfactor)**np.arange(pad) +ypad = (np.ones(pad)*padfactor)**np.arange(pad) + +core = 15 +xcore = np.ones(core) +ycore = np.ones(core) + +h1 = np.r_[xpad[::-1],xcore,xpad] +h2 = np.r_[ypad[::-1],ycore,ypad] + +mesh = TensorMesh([h1, h2]) +mesh.plotGrid() +plt.axis('tight') + +plt.show() diff --git a/docs/examples/mesh/plot_grid_2D.py b/docs/examples/mesh/plot_grid_2D.py index c53a8e02..258fbdfc 100644 --- a/docs/examples/mesh/plot_grid_2D.py +++ b/docs/examples/mesh/plot_grid_2D.py @@ -5,7 +5,7 @@ from SimPEG import TensorMesh h1 = np.linspace(.1,.5,3) h2 = np.linspace(.1,.5,5) mesh = TensorMesh([h1, h2]) -mesh.plotGrid() +mesh.plotGrid(nodes=True, faces=True, centers=True, lines=True) plt.show() diff --git a/docs/examples/mesh/plot_grid_3D.py b/docs/examples/mesh/plot_grid_3D.py index 8da62346..caf8dbd2 100644 --- a/docs/examples/mesh/plot_grid_3D.py +++ b/docs/examples/mesh/plot_grid_3D.py @@ -6,7 +6,7 @@ h1 = np.linspace(.1,.5,3) h2 = np.linspace(.1,.5,5) h3 = np.linspace(.1,.5,3) mesh = TensorMesh([h1,h2,h3]) -mesh.plotGrid() +mesh.plotGrid(nodes=True, faces=True, centers=True, lines=True) plt.show() From 71668deb1b5f8ae50eb8060c842986ec1b44d9dd Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 30 Sep 2013 16:46:02 -0700 Subject: [PATCH 3/4] inner product documentation --- SimPEG/InnerProducts.py | 197 +++++++++++++++++++++++++++++- SimPEG/LogicallyOrthogonalMesh.py | 17 ++- docs/api_LOMView.rst | 4 +- docs/index.rst | 2 +- 4 files changed, 204 insertions(+), 16 deletions(-) diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py index 1add3e72..18803344 100644 --- a/SimPEG/InnerProducts.py +++ b/SimPEG/InnerProducts.py @@ -43,6 +43,88 @@ class InnerProducts(object): def getFaceInnerProduct(mesh, mu=None, returnP=False): + """ + :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :rtype: scipy.csr_matrix + :return: M, the inner product matrix + + Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows: + + .. math:: + \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{1} & 0 \\\\ 0 & 0 & \mu_{1} \end{matrix}\\right] + + \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{2} & 0 \\\\ 0 & 0 & \mu_{3} \end{matrix}\\right] + + \left[\\begin{matrix} \mu_{1} & \mu_{4} & \mu_{5} \\\\ \mu_{4} & \mu_{2} & \mu_{6} \\\\ \mu_{5} & \mu_{6} & \mu_{3} \end{matrix}\\right] + + Example problem for DC resistivity: + + .. math:: + + \sigma^{-1}\mathbf{J} = \\nabla \phi + + We can define in weak form by integrating with a general face function F: + + .. math:: + + \int_{\\text{cell}}{\sigma^{-1}\mathbf{J} \cdot \mathbf{F}} = \int_{\\text{cell}}{\\nabla \phi \cdot \mathbf{F}} + + \int_{\\text{cell}}{\sigma^{-1}\mathbf{J} \cdot \mathbf{F}} = \int_{\\text{cell}}{(\\nabla \cdot \mathbf{F}) \phi } + \int_{\partial \\text{cell}}{ \phi \mathbf{F} \cdot \mathbf{n}} + + We can then discretize for every cell: + + .. math:: + + v_{\\text{cell}} \sigma^{-1} (\mathbf{J}_x \mathbf{F}_x +\mathbf{J}_y \mathbf{F}_y + \mathbf{J}_z \mathbf{F}_z ) = -\phi^{\\top} v_{\\text{cell}} (\mathbf{D}_{\\text{cell}} \mathbf{F}) + \\text{BC} + + We can represent this in vector form (again this is for every cell), and will generalize for the case of anisotropic (tensor) sigma. + + .. math:: + + \mathbf{F}_c^{\\top} (\sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}}) \mathbf{J}_c = -\phi^{\\top} v_{\\text{cell}}( v_\\text{cell}^{-1} \mathbf{D}_{\\text{cell}} \mathbf{A} \mathbf{F}) + \\text{BC} + + We multiply by volume on each side of the tensor conductivity to keep symmetry in the system. Here J_c is the Cartesian J (on the faces) and must be calculated differently depending on the mesh: + + .. math:: + \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\\text{TENSOR} = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\\text{LOM} + + Here the i index refers to where we choose to approximate this integral. We will approximate this relation at every node of the cell, there are 8 in 3D, using a projection matrix Q_i to pick the appropriate fluxes. We will then average to the cell center: + + .. math:: + + \mathbf{F}^{\\top} + {1\over 8} + \left(\sum_{i=1}^8 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + \mathbf{J} + = + -\mathbf{F}^{\\top} \mathbf{A} \mathbf{D}_{\\text{cell}}^{\\top} \phi + \\text{BC} + + \mathbf{M}(\Sigma^{-1}) \mathbf{J} + = + -\mathbf{A} \mathbf{D}_{\\text{cell}}^{\\top} \phi + \\text{BC} + + \mathbf{M}(\Sigma^{-1}) = {1\over 8} + \left(\sum_{i=1}^8 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + The M is returned if mu is set equal to \Sigma^{-1}. + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P000, P001, P010, P011, P100, P101, P110, P111] + + Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + + """ if mu is None: # default is ones mu = np.ones((mesh.nC, 1)) @@ -82,10 +164,10 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False): # 100 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k # 010 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k # 110 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k - # 001 | i ,j ,k | i , j, k | i, j , k | i, j, k+1 - # 101 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k+1 - # 011 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k+1 - # 111 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k+1 + # 001 | i ,j ,k+1 | i , j, k | i, j , k | i, j, k+1 + # 101 | i+1,j ,k+1 | i+1, j, k | i, j , k | i, j, k+1 + # 011 | i ,j+1,k+1 | i , j, k | i, j+1, k | i, j, k+1 + # 111 | i+1,j+1,k+1 | i+1, j, k | i, j+1, k | i, j, k+1 # Square root of cell volume multiplied by 1/8 v = np.sqrt(0.125*mesh.vol) @@ -120,6 +202,42 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False): def getFaceInnerProduct2D(mesh, mu=None, returnP=False): + """ + :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 2, or 3)) + :param bool returnP: returns the projection matrices + :rtype: scipy.csr_matrix + :return: M, the inner product matrix + + Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows: + + .. math:: + \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{1} \end{matrix}\\right] + + \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{2} \end{matrix}\\right] + + \left[\\begin{matrix} \mu_{1} & \mu_{3} \\\\ \mu_{3} & \mu_{2} \end{matrix}\\right] + + + .. math:: + + \mathbf{M}(\Sigma^{-1}) = {1\over 4} + \left(\sum_{i=1}^4 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P00, P10, P01, P11] + + Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + + """ if mu is None: # default is ones mu = np.ones((mesh.nC, 1)) @@ -185,6 +303,41 @@ def getFaceInnerProduct2D(mesh, mu=None, returnP=False): def getEdgeInnerProduct(mesh, sigma=None, returnP=False): + """ + :param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :rtype: scipy.csr_matrix + :return: M, the inner product matrix + + + Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows: + + .. math:: + \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{1} & 0 \\\\ 0 & 0 & \sigma_{1} \end{matrix}\\right] + + \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{2} & 0 \\\\ 0 & 0 & \sigma_{3} \end{matrix}\\right] + + \left[\\begin{matrix} \sigma_{1} & \sigma_{4} & \sigma_{5} \\\\ \sigma_{4} & \sigma_{2} & \sigma_{6} \\\\ \sigma_{5} & \sigma_{6} & \sigma_{3} \end{matrix}\\right] + + What is returned: + + .. math:: + \mathbf{M}(\Sigma^{-1}) = {1\over 8} + \left(\sum_{i=1}^8 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P000, P001, P010, P011, P100, P101, P110, P111] + + Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + """ if sigma is None: # default is ones sigma = np.ones((mesh.nC, 1)) @@ -262,6 +415,42 @@ def getEdgeInnerProduct(mesh, sigma=None, returnP=False): def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False): + """ + :param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 2, or 3)) + :param bool returnP: returns the projection matrices + :rtype: scipy.csr_matrix + :return: M, the inner product matrix + + Depending on the number of columns (either 1, 2, or 3) of sigma, the material property is interpreted as follows: + + .. math:: + \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{1} \end{matrix}\\right] + + \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{2} \end{matrix}\\right] + + \left[\\begin{matrix} \sigma_{1} & \sigma_{3} \\\\ \sigma_{3} & \sigma_{2} \end{matrix}\\right] + + + .. math:: + + \mathbf{M}(\Sigma^{-1}) = {1\over 4} + \left(\sum_{i=1}^4 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P00, P10, P01, P11] + + Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + + """ if sigma is None: # default is ones sigma = np.ones((mesh.nC, 1)) diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/LogicallyOrthogonalMesh.py index 49e2ee07..b9516a6d 100644 --- a/SimPEG/LogicallyOrthogonalMesh.py +++ b/SimPEG/LogicallyOrthogonalMesh.py @@ -258,17 +258,16 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): area = property(**area()) def normals(): - doc = """ -Face normals: calling this will average -the computed normals so that there is one -per face. This is especially relevant in -3D, as there are up to 4 different normals -for each face that will be different. + doc = """Face normals: calling this will average + the computed normals so that there is one + per face. This is especially relevant in + 3D, as there are up to 4 different normals + for each face that will be different. -To reshape the normals into a matrix and get the y component: + To reshape the normals into a matrix and get the y component: -NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M') -""" + NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M') + """ def fget(self): if(self._normals is None): diff --git a/docs/api_LOMView.rst b/docs/api_LOMView.rst index 67e1bb91..61630c26 100644 --- a/docs/api_LOMView.rst +++ b/docs/api_LOMView.rst @@ -1,8 +1,8 @@ .. _api_LOMView: LOM View -*********** +******** -.. automodule:: SimPEG.LOMView +.. automodule:: SimPEG.LomView :members: :undoc-members: diff --git a/docs/index.rst b/docs/index.rst index ce50bc3d..22c9d5e9 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -26,7 +26,7 @@ Meshing & Operators api_TensorMesh api_TensorView api_LogicallyOrthogonalMesh - api_LomView + api_LOMView api_DiffOperators api_InnerProducts From 063a3d8b57b7f18462f2266cbec64bb38150d2ae Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 30 Sep 2013 17:19:00 -0700 Subject: [PATCH 4/4] updates to documentation --- SimPEG/DiffOperators.py | 12 +-- SimPEG/InnerProducts.py | 182 +++++++++++++++++++++++----------------- SimPEG/LomView.py | 10 ++- 3 files changed, 119 insertions(+), 85 deletions(-) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index f8c1feb9..9a8a9db9 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -111,12 +111,12 @@ class DiffOperators(object): """ Function that sets the boundary conditions for cell-centred derivative operators. - Examples: + 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 + 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): @@ -266,7 +266,7 @@ class DiffOperators(object): nodalAve = property(**nodalAve()) def nodalVectorAve(): - doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension seperate." + doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension separate." def fget(self): if(self._nodalVectorAve is None): diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py index 18803344..fca632f4 100644 --- a/SimPEG/InnerProducts.py +++ b/SimPEG/InnerProducts.py @@ -6,59 +6,11 @@ import numpy as np class InnerProducts(object): """ Class creates the inner product matrices that you need! - """ - def __init__(self): - raise Exception('InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.') - def getFaceInnerProduct(self, mu=None, returnP=False): - if self.dim == 2: - return getFaceInnerProduct2D(self, mu, returnP) - elif self.dim == 3: - return getFaceInnerProduct(self, mu, returnP) - - def getEdgeInnerProduct(self, sigma=None, returnP=False): - if self.dim == 2: - return getEdgeInnerProduct2D(self, sigma, returnP) - elif self.dim == 3: - return getEdgeInnerProduct(self, sigma, returnP) - -# ------------------------ Geometries ------------------------------ -# -# -# node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1) -# / / -# / / | -# edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k) -# / / | -# / / | -# node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k) -# | | | -# | | node(i+1,j+1,k+1) -# | | / -# edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k) -# | | / -# | | / -# | |/ -# node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) + InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class. -def getFaceInnerProduct(mesh, mu=None, returnP=False): - """ - :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param bool returnP: returns the projection matrices - :rtype: scipy.csr_matrix - :return: M, the inner product matrix - - Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows: - - .. math:: - \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{1} & 0 \\\\ 0 & 0 & \mu_{1} \end{matrix}\\right] - - \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{2} & 0 \\\\ 0 & 0 & \mu_{3} \end{matrix}\\right] - - \left[\\begin{matrix} \mu_{1} & \mu_{4} & \mu_{5} \\\\ \mu_{4} & \mu_{2} & \mu_{6} \\\\ \mu_{5} & \mu_{6} & \mu_{3} \end{matrix}\\right] - - Example problem for DC resistivity: + **Example problem for DC resistivity** .. math:: @@ -89,14 +41,16 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False): .. math:: \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\\text{TENSOR} = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\\text{LOM} - Here the i index refers to where we choose to approximate this integral. We will approximate this relation at every node of the cell, there are 8 in 3D, using a projection matrix Q_i to pick the appropriate fluxes. We will then average to the cell center: + Here the i index refers to where we choose to approximate this integral. + We will approximate this relation at every node of the cell, there are 8 in 3D, using a projection matrix Q_i to pick the appropriate fluxes. + We will then average to the cell center. For the TENSOR mesh, this looks like: .. math:: \mathbf{F}^{\\top} {1\over 8} \left(\sum_{i=1}^8 - \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \mathbf{Q}_{(i)}^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{Q}_{(i)} \\right) \mathbf{J} = @@ -108,16 +62,92 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False): \mathbf{M}(\Sigma^{-1}) = {1\over 8} \left(\sum_{i=1}^8 - \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \mathbf{Q}_{(i)}^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{Q}_{(i)} \\right) The M is returned if mu is set equal to \Sigma^{-1}. + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes). + Here each P (3*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + """ + def __init__(self): + raise Exception('InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.') + + def getFaceInnerProduct(self, mu=None, returnP=False): + """Wrapper function, + + :py:func:`SimPEG.InnerProducts.getEdgeInnerProduct` + + :py:func:`SimPEG.InnerProducts.getEdgeInnerProduct2D` + """ + if self.dim == 2: + return getFaceInnerProduct2D(self, mu, returnP) + elif self.dim == 3: + return getFaceInnerProduct(self, mu, returnP) + + def getEdgeInnerProduct(self, sigma=None, returnP=False): + """Wrapper function, + + :py:func:`SimPEG.InnerProducts.getFaceInnerProduct` + + :py:func:`SimPEG.InnerProducts.getFaceInnerProduct2D` + """ + if self.dim == 2: + return getEdgeInnerProduct2D(self, sigma, returnP) + elif self.dim == 3: + return getEdgeInnerProduct(self, sigma, returnP) + +# ------------------------ Geometries ------------------------------ +# +# +# node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1) +# / / +# / / | +# edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k) +# / / | +# / / | +# node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k) +# | | | +# | | node(i+1,j+1,k+1) +# | | / +# edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k) +# | | / +# | | / +# | |/ +# node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) + + +def getFaceInnerProduct(mesh, mu=None, returnP=False): + """ + :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (sum(nF), sum(nF)) + + Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows: + + .. math:: + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{1} & 0 \\\\ 0 & 0 & \mu_{1} \end{matrix}\\right] + + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{2} & 0 \\\\ 0 & 0 & \mu_{3} \end{matrix}\\right] + + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{4} & \mu_{5} \\\\ \mu_{4} & \mu_{2} & \mu_{6} \\\\ \mu_{5} & \mu_{6} & \mu_{3} \end{matrix}\\right] + + \mathbf{M}(\\vec{\mu}) = {1\over 8} + \left(\sum_{i=1}^8 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: P = [P000, P001, P010, P011, P100, P101, P110, P111] - Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + Here each P (3*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: .. math:: \mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} @@ -206,23 +236,23 @@ def getFaceInnerProduct2D(mesh, mu=None, returnP=False): :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 2, or 3)) :param bool returnP: returns the projection matrices :rtype: scipy.csr_matrix - :return: M, the inner product matrix + :return: M, the inner product matrix (sum(nF), sum(nF)) Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows: .. math:: - \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{1} \end{matrix}\\right] + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{1} \end{matrix}\\right] - \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{2} \end{matrix}\\right] + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{2} \end{matrix}\\right] - \left[\\begin{matrix} \mu_{1} & \mu_{3} \\\\ \mu_{3} & \mu_{2} \end{matrix}\\right] + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{3} \\\\ \mu_{3} & \mu_{2} \end{matrix}\\right] .. math:: - \mathbf{M}(\Sigma^{-1}) = {1\over 4} + \mathbf{M}(\\vec{\mu}) = {1\over 4} \left(\sum_{i=1}^4 - \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}} \mathbf{J}_c \\right) @@ -230,7 +260,7 @@ def getFaceInnerProduct2D(mesh, mu=None, returnP=False): P = [P00, P10, P01, P11] - Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + Here each P (2*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: .. math:: \mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} @@ -307,31 +337,31 @@ def getEdgeInnerProduct(mesh, sigma=None, returnP=False): :param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :param bool returnP: returns the projection matrices :rtype: scipy.csr_matrix - :return: M, the inner product matrix + :return: M, the inner product matrix (sum(nE), sum(nE)) - Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows: + Depending on the number of columns (either 1, 3, or 6) of sigma, the material property is interpreted as follows: .. math:: - \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{1} & 0 \\\\ 0 & 0 & \sigma_{1} \end{matrix}\\right] + \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{1} & 0 \\\\ 0 & 0 & \sigma_{1} \end{matrix}\\right] - \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{2} & 0 \\\\ 0 & 0 & \sigma_{3} \end{matrix}\\right] + \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{2} & 0 \\\\ 0 & 0 & \sigma_{3} \end{matrix}\\right] - \left[\\begin{matrix} \sigma_{1} & \sigma_{4} & \sigma_{5} \\\\ \sigma_{4} & \sigma_{2} & \sigma_{6} \\\\ \sigma_{5} & \sigma_{6} & \sigma_{3} \end{matrix}\\right] + \Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{4} & \sigma_{5} \\\\ \sigma_{4} & \sigma_{2} & \sigma_{6} \\\\ \sigma_{5} & \sigma_{6} & \sigma_{3} \end{matrix}\\right] What is returned: .. math:: - \mathbf{M}(\Sigma^{-1}) = {1\over 8} + \mathbf{M}(\Sigma) = {1\over 8} \left(\sum_{i=1}^8 - \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}} \mathbf{J}_c \\right) If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: P = [P000, P001, P010, P011, P100, P101, P110, P111] - Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + Here each P (3*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: .. math:: \mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} @@ -419,23 +449,23 @@ def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False): :param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 2, or 3)) :param bool returnP: returns the projection matrices :rtype: scipy.csr_matrix - :return: M, the inner product matrix + :return: M, the inner product matrix (sum(nE), sum(nE)) Depending on the number of columns (either 1, 2, or 3) of sigma, the material property is interpreted as follows: .. math:: - \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{1} \end{matrix}\\right] + \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{1} \end{matrix}\\right] - \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{2} \end{matrix}\\right] + \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{2} \end{matrix}\\right] - \left[\\begin{matrix} \sigma_{1} & \sigma_{3} \\\\ \sigma_{3} & \sigma_{2} \end{matrix}\\right] + \Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{3} \\\\ \sigma_{3} & \sigma_{2} \end{matrix}\\right] .. math:: - \mathbf{M}(\Sigma^{-1}) = {1\over 4} + \mathbf{M}(\Sigma) = {1\over 4} \left(\sum_{i=1}^4 - \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}} \mathbf{J}_c \\right) @@ -443,7 +473,7 @@ def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False): P = [P00, P10, P01, P11] - Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + Here each P (2*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: .. math:: \mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} diff --git a/SimPEG/LomView.py b/SimPEG/LomView.py index 4c9b1dab..4b4f36a3 100644 --- a/SimPEG/LomView.py +++ b/SimPEG/LomView.py @@ -7,15 +7,19 @@ from utils import mkvc class LomView(object): """ - Provides viewing functions for TensorMesh + Provides viewing functions for LogicallyOrthogonalMesh + + This class is inherited by LogicallyOrthogonalMesh - This class is inherited by TensorMesh """ def __init__(self): pass def plotGrid(self, length=0.05): - """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.""" + """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. + + .. plot:: examples/mesh/plot_LogicallyOrthogonalMesh.py + """ NN = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: fig = plt.figure(2)