From 1c17c3341e54d6ca1addd713998b1dd51176f884 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 26 Feb 2014 15:12:30 -0800 Subject: [PATCH 01/24] Major updates to innerProducts documentation --- SimPEG/Mesh/InnerProducts.py | 295 ++++++----------------------------- docs/api_InnerProducts.rst | 164 +++++++++++++++++++ docs/api_Mesh.rst | 8 - docs/index.rst | 1 + 4 files changed, 215 insertions(+), 253 deletions(-) create mode 100644 docs/api_InnerProducts.rst diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 28befd6e..87d106d2 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -5,183 +5,51 @@ import numpy as np class InnerProducts(object): """ - Class creates the inner product matrices that you need! - - InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class. - - - **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. For the TENSOR mesh, this looks like: - - .. math:: - - \mathbf{F}^{\\top} - {1\over 8} - \left(\sum_{i=1}^8 - \mathbf{Q}_{(i)}^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{Q}_{(i)} - \\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{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. + This is a base for the SimPEG.Mesh classes. This mixIn creates the all 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(M, mu=None, returnP=False): + def getFaceInnerProduct(self, materialProperty=None, returnP=False): """ - :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :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 :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, P100, P010, P110, P001, P101, P011, P111] - - 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. - - **For 2D:** - - Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows: - - .. math:: - \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{1} \end{matrix}\\right] - - \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{2} \end{matrix}\\right] - - \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{3} \\\\ \mu_{3} & \mu_{2} \end{matrix}\\right] - - - .. math:: - - \mathbf{M}(\\vec{\mu}) = {1\over 4} - \left(\sum_{i=1}^4 - \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 = [P00, P10, P01, P11] - - 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)} - - Note that this is completed for each cell in the mesh at the same time. - + :return: M, the inner product matrix (nF, nF) """ - if M.dim == 1: - v = np.sqrt(0.5*M.vol) - V1 = sdiag(v) # We will multiply on each side to keep symmetry + d = self.dim + # We will multiply by sqrt on each side to keep symmetry + V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol))) - Px = _getFacePx(M) - P000 = V1*Px('fXm') - P100 = V1*Px('fXp') - elif M.dim == 2: - # Square root of cell volume multiplied by 1/4 - v = np.sqrt(0.25*M.vol) - V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry + if d == 1: + Px = _getFacePx(self) + P000 = V*Px('fXm') + P100 = V*Px('fXp') + elif d == 2: + Pxx = _getFacePxx(self) + P000 = V*Pxx('fXm', 'fYm') + P100 = V*Pxx('fXp', 'fYm') + P010 = V*Pxx('fXm', 'fYp') + P110 = V*Pxx('fXp', 'fYp') + elif d == 3: + Pxxx = _getFacePxxx(self) + P000 = V*Pxxx('fXm', 'fYm', 'fZm') + P100 = V*Pxxx('fXp', 'fYm', 'fZm') + P010 = V*Pxxx('fXm', 'fYp', 'fZm') + P110 = V*Pxxx('fXp', 'fYp', 'fZm') + P001 = V*Pxxx('fXm', 'fYm', 'fZp') + P101 = V*Pxxx('fXp', 'fYm', 'fZp') + P011 = V*Pxxx('fXm', 'fYp', 'fZp') + P111 = V*Pxxx('fXp', 'fYp', 'fZp') - Pxx = _getFacePxx(M) - P000 = V2*Pxx('fXm', 'fYm') - P100 = V2*Pxx('fXp', 'fYm') - P010 = V2*Pxx('fXm', 'fYp') - P110 = V2*Pxx('fXp', 'fYp') - elif M.dim == 3: - # Square root of cell volume multiplied by 1/8 - v = np.sqrt(0.125*M.vol) - V3 = sdiag(np.r_[v, v, v]) # We will multiply on each side to keep symmetry - - Pxxx = _getFacePxxx(M) - P000 = V3*Pxxx('fXm', 'fYm', 'fZm') - P100 = V3*Pxxx('fXp', 'fYm', 'fZm') - P010 = V3*Pxxx('fXm', 'fYp', 'fZm') - P110 = V3*Pxxx('fXp', 'fYp', 'fZm') - P001 = V3*Pxxx('fXm', 'fYm', 'fZp') - P101 = V3*Pxxx('fXp', 'fYm', 'fZp') - P011 = V3*Pxxx('fXm', 'fYp', 'fZp') - P111 = V3*Pxxx('fXp', 'fYp', 'fZp') - - Mu = makePropertyTensor(M, mu) + Mu = makePropertyTensor(self, materialProperty) A = P000.T*Mu*P000 + P100.T*Mu*P100 P = [P000, P100] - if M.dim > 1: + if d > 1: A = A + P010.T*Mu*P010 + P110.T*Mu*P110 P += [P010, P110] - if M.dim > 2: + if d > 2: A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111 P += [P001, P101, P011, P111] if returnP: @@ -189,91 +57,28 @@ class InnerProducts(object): else: return A - def getEdgeInnerProduct(M, sigma=None, returnP=False): + def getEdgeInnerProduct(self, materialProperty=None, returnP=False): """ - :param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :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 :rtype: scipy.csr_matrix - :return: M, the inner product matrix (sum(nE), sum(nE)) - - - Depending on the number of columns (either 1, 3, or 6) of sigma, the material property is interpreted as follows: - - .. math:: - \Sigma = \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_{2} & 0 \\\\ 0 & 0 & \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\over 8} - \left(\sum_{i=1}^8 - \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, P100, P010, P110, P001, P101, P011, P111] - - 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)} - - Note that this is completed for each cell in the mesh at the same time. - - **For 2D:** - - Depending on the number of columns (either 1, 2, or 3) of sigma, the material property is interpreted as follows: - - .. math:: - \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{1} \end{matrix}\\right] - - \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \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\over 4} - \left(\sum_{i=1}^4 - \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 = [P00, P10, P01, P11] - - 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)} - - Note that this is completed for each cell in the mesh at the same time. - + :return: M, the inner product matrix (nE, nE) """ - if M.dim == 1: + d = self.dim + # We will multiply by sqrt on each side to keep symmetry + V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol))) + + if d == 1: raise NotImplementedError('getEdgeInnerProduct not implemented for 1D') # We will multiply by V on each side to keep symmetry - elif M.dim == 2: - # Square root of cell volume multiplied by 1/4 - v = np.sqrt(0.25*M.vol) - V = sdiag(np.r_[v, v]) - eP = _getEdgePxx(M) + elif d == 2: + eP = _getEdgePxx(self) P000 = V*eP('eX0', 'eY0') P100 = V*eP('eX0', 'eY1') P010 = V*eP('eX1', 'eY0') P110 = V*eP('eX1', 'eY1') - elif M.dim == 3: - # Square root of cell volume multiplied by 1/8 - v = np.sqrt(0.125*M.vol) - V = sdiag(np.r_[v, v, v]) - eP = _getEdgePxxx(M) + elif d == 3: + eP = _getEdgePxxx(self) P000 = V*eP('eX0', 'eY0', 'eZ0') P100 = V*eP('eX0', 'eY1', 'eZ1') P010 = V*eP('eX1', 'eY0', 'eZ2') @@ -283,11 +88,11 @@ class InnerProducts(object): P011 = V*eP('eX3', 'eY2', 'eZ2') P111 = V*eP('eX3', 'eY3', 'eZ3') - Sigma = makePropertyTensor(M, sigma) - A = P000.T*Sigma*P000 + P100.T*Sigma*P100 + P010.T*Sigma*P010 + P110.T*Sigma*P110 + Mu = makePropertyTensor(self, materialProperty) + A = P000.T*Mu*P000 + P100.T*Mu*P100 + P010.T*Mu*P010 + P110.T*Mu*P110 P = [P000, P100, P010, P110] - if M.dim == 3: - A = A + P001.T*Sigma*P001 + P101.T*Sigma*P101 + P011.T*Sigma*P011 + P111.T*Sigma*P111 + if d == 3: + A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111 P += [P001, P101, P011, P111] if returnP: return A, P @@ -381,11 +186,11 @@ def _getFacePxx_Rectangular(M): 0 1 f2(Ym) - Pxx('m','m') = | 1, 0, 0, 0 | - | 0, 0, 1, 0 | + Pxx('fXm','fYm') = | 1, 0, 0, 0 | + | 0, 0, 1, 0 | - Pxx('p','m') = | 0, 1, 0, 0 | - | 0, 0, 1, 0 | + Pxx('fXp','fYm') = | 0, 1, 0, 0 | + | 0, 0, 1, 0 | """ i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy)) diff --git a/docs/api_InnerProducts.rst b/docs/api_InnerProducts.rst new file mode 100644 index 00000000..0909f086 --- /dev/null +++ b/docs/api_InnerProducts.rst @@ -0,0 +1,164 @@ +.. _api_InnerProducts: + + +Inner Products +************** + +By using the weak formulation of many of the PDEs in geophysical applications, we can rapidly develop discretizations. Much of this work, however, needs a good understanding of how to approximate inner products on our discretized meshes. We will define the inner product as: + +.. math:: + + \left(a,b\right) = \int_\Omega{a \cdot b}{\partial v} + +where a and b are either scalars or vectors. + +.. note:: + + The InnerProducts class is a base class providing inner product matrices for meshes and cannot run on its own. + + +Example problem for DC resistivity +---------------------------------- +We will start with the formulation of the Direct Current (DC) resistivity problem in geophysics. + + +.. math:: + + \frac{1}{\sigma}\vec{j} = \nabla \phi \\ + + \nabla\cdot \vec{j} = q + +In the following discretization, \\\( \\sigma \\\) and \\\( \\phi \\\) +will be discretized on the cell-centers and the flux, \\\(\\vec{j}\\\), +will be on the faces. We will use the weak formulation to discretize +the DC resistivity equation. + +We can define in weak form by integrating with a general face function \\\(\\vec{f}\\\): + +.. math:: + + \int_{\Omega}{\sigma^{-1}\vec{j} \cdot \vec{f}} = \int_{\Omega}{\nabla \phi \cdot \vec{f}} + +Here we can integrate the right side by parts, + +.. math:: + + \nabla\cdot(\phi\vec{f})=\nabla\phi\cdot\vec{f} + \phi\nabla\cdot\vec{f} + +and rearrange it, and apply the Divergence theorem. + +.. math:: + + \int_{\Omega}{\sigma^{-1}\vec{j} \cdot \vec{f}} = + - \int_{\Omega}{(\phi \nabla \cdot \vec{f})} + + \int_{\partial \Omega}{ \phi \vec{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} + +.. note:: + + We have discretized the dot product above, but remember that we do not really have a single vector \\\(\\mathbf{J}\\\), but approximations of \\\(\\vec{j}\\\) on each face of our cell. In 2D that means 2 approximations of \\\(\\mathbf{J}_x\\\) and 2 approximations of \\\(\\mathbf{J}_y\\\). In 3D we also have 2 approximations of \\\(\\mathbf{J}_z\\\). + +Regardless of how we choose to approximate this dot product, 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 square-root of volume on each side of the tensor conductivity to keep symmetry in the system. Here \\\(\\mathbf{J}_c\\\) is the Cartesian \\\(\\mathbf{J}\\\) (on the faces that we choose to use in our approximation) and must be calculated differently depending on the mesh: + +.. math:: + \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\text{TENSOR} \\ + \mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{LRM} + +Here the \\\(i\\\) index refers to where we choose to approximate this integral, as discussed in the note above. +We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix \\\( \\mathbf{Q}_{(i)} \\\) to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like: + +.. math:: + + \mathbf{F}^{\top} + {1\over 8} + \left(\sum_{i=1}^8 + \mathbf{Q}_{(i)}^{\top} \sqrt{v_{\text{cell}}} \Sigma^{-1} \sqrt{v_{\text{cell}}} \mathbf{Q}_{(i)} + \right) + \mathbf{J} + = + -\mathbf{F}^{\top} \mathbf{A} \mathbf{D}_{\text{cell}}^{\top} \phi + \text{BC} + +Or, when generalizing to the entire mesh and dropping our general face function: + +.. math:: + + \mathbf{M}^f_{\Sigma^{-1}} \mathbf{J} + = + -\mathbf{A} \mathbf{D}^{\top} \phi + \text{BC} + +By defining the faceInnerProduct in 3D (8 combinations of fluxes) to be: + +.. math:: + + \mathbf{M}^f_{\Sigma^{-1}} = {1\over 8} + \left(\sum_{i=1}^8 + \mathbf{P}_{(i)}^{\top} \Sigma^{-1} \mathbf{P}_{(i)} + \right) + +The M is returned when given the input of \\\( \\Sigma^{-1} \\\). + +Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of the projection, volume, and any normalization to Cartesian coordinates (where the dot product is well defined): + +.. math:: + + \mathbf{P}_{(i)} = \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{LRM only}} \mathbf{Q}_{(i)} + +Where \\\(d\\\) is the dimension of the mesh. + +.. note:: + + This is actually completed for each cell in the mesh at the same time, and the full matrices are returned. + +If ``returnP=True`` is requested in any of these methods the projection matrices are returned as a list ordered by nodes around which the fluxes were approximated:: + + # In 3D + P = [P000, P100, P010, P110, P001, P101, P011, P111] + # In 2D + P = [P00, P10, P01, P11] + +The derivation for ``edgeInnerProducts`` is exactly the same, however, when we approximate the integral using the fields around each node, the projection matrices look a bit different because we have 12 edges in 3D instead of just 6 faces. The interface to the code is exactly the same. + +Defining Tensor Properties +-------------------------- + +**For 3D:** + +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] + +**For 2D:** + + Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows: + +.. math:: + \vec{\mu} = \left[\begin{matrix} \mu_{1} & 0 \\ 0 & \mu_{1} \end{matrix}\right] + + \vec{\mu} = \left[\begin{matrix} \mu_{1} & 0 \\ 0 & \mu_{2} \end{matrix}\right] + + \vec{\mu} = \left[\begin{matrix} \mu_{1} & \mu_{3} \\ \mu_{3} & \mu_{2} \end{matrix}\right] + + +The API +------- + +.. automodule:: SimPEG.Mesh.InnerProducts + :members: + :undoc-members: diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst index a964b24d..df5f31d8 100644 --- a/docs/api_Mesh.rst +++ b/docs/api_Mesh.rst @@ -29,7 +29,6 @@ Logically Orthogonal Mesh :inherited-members: - Base Mesh ========= @@ -38,13 +37,6 @@ Base Mesh :undoc-members: -Inner Products -============== - -.. automodule:: SimPEG.Mesh.InnerProducts - :members: - :undoc-members: - Differential Operators ====================== diff --git a/docs/index.rst b/docs/index.rst index 84e0f925..6a282e36 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -30,6 +30,7 @@ Meshing & Operators :maxdepth: 2 api_Mesh + api_InnerProducts Forward Problems **************** From 0dec95aa4ffeaaea4383a23c82fc0c30f056ae30 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 26 Feb 2014 15:38:28 -0800 Subject: [PATCH 02/24] updates/typos in docs --- docs/api_InnerProducts.rst | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/docs/api_InnerProducts.rst b/docs/api_InnerProducts.rst index 0909f086..dec9c4c8 100644 --- a/docs/api_InnerProducts.rst +++ b/docs/api_InnerProducts.rst @@ -57,7 +57,7 @@ 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} + 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} .. note:: @@ -67,7 +67,9 @@ Regardless of how we choose to approximate this dot product, we can represent th .. 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} + \mathbf{F}_c^{\top} (\sqrt{v_{\text{cell}}} \Sigma^{-1} \sqrt{v_{\text{cell}}}) \mathbf{J}_c = + -\phi^{\top} v_{\text{cell}} \mathbf{D}_{\text{cell}} \mathbf{F}) + + \text{BC} We multiply by square-root of volume on each side of the tensor conductivity to keep symmetry in the system. Here \\\(\\mathbf{J}_c\\\) is the Cartesian \\\(\\mathbf{J}\\\) (on the faces that we choose to use in our approximation) and must be calculated differently depending on the mesh: @@ -87,7 +89,7 @@ We will approximate this integral by taking the fluxes clustered around every no \right) \mathbf{J} = - -\mathbf{F}^{\top} \mathbf{A} \mathbf{D}_{\text{cell}}^{\top} \phi + \text{BC} + -\mathbf{F}^{\top} \mathbf{D}_{\text{cell}}^{\top} v_{\text{cell}} \phi + \text{BC} Or, when generalizing to the entire mesh and dropping our general face function: @@ -95,18 +97,18 @@ Or, when generalizing to the entire mesh and dropping our general face function: \mathbf{M}^f_{\Sigma^{-1}} \mathbf{J} = - -\mathbf{A} \mathbf{D}^{\top} \phi + \text{BC} + - \mathbf{D}^{\top} \text{diag}(\mathbf{v}) \phi + \text{BC} -By defining the faceInnerProduct in 3D (8 combinations of fluxes) to be: +By defining the faceInnerProduct (8 combinations of fluxes in 3D, 4 in 2D, 2 in 1D) to be: .. math:: - \mathbf{M}^f_{\Sigma^{-1}} = {1\over 8} - \left(\sum_{i=1}^8 + \mathbf{M}^f_{\Sigma^{-1}} = + \sum_{i=1}^{2^d} \mathbf{P}_{(i)}^{\top} \Sigma^{-1} \mathbf{P}_{(i)} - \right) -The M is returned when given the input of \\\( \\Sigma^{-1} \\\). +Where \\\(d\\\) is the dimension of the mesh. +The \\\( \\mathbf{M}^f \\\) is returned when given the input of \\\( \\Sigma^{-1} \\\). Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of the projection, volume, and any normalization to Cartesian coordinates (where the dot product is well defined): @@ -114,8 +116,6 @@ Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of \mathbf{P}_{(i)} = \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{LRM only}} \mathbf{Q}_{(i)} -Where \\\(d\\\) is the dimension of the mesh. - .. note:: This is actually completed for each cell in the mesh at the same time, and the full matrices are returned. @@ -126,6 +126,8 @@ If ``returnP=True`` is requested in any of these methods the projection matrices P = [P000, P100, P010, P110, P001, P101, P011, P111] # In 2D P = [P00, P10, P01, P11] + # In 1D + P = [P0, P1] The derivation for ``edgeInnerProducts`` is exactly the same, however, when we approximate the integral using the fields around each node, the projection matrices look a bit different because we have 12 edges in 3D instead of just 6 faces. The interface to the code is exactly the same. From b0b4a00172d082172f975c40aa15b0b3b88d9838 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 26 Feb 2014 16:25:40 -0800 Subject: [PATCH 03/24] docs updates and notation changes --- SimPEG/Mesh/InnerProducts.py | 35 +++++++++++++++++------------------ docs/api_DiffOps.rst | 8 ++++++++ docs/api_Mesh.rst | 9 +-------- docs/api_Solver.rst | 10 ++++++++++ docs/api_Utils.rst | 8 -------- docs/index.rst | 36 +++++++++++++++--------------------- 6 files changed, 51 insertions(+), 55 deletions(-) create mode 100644 docs/api_DiffOps.rst create mode 100644 docs/api_Solver.rst diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 87d106d2..a797a4a9 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -22,25 +22,25 @@ class InnerProducts(object): V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol))) if d == 1: - Px = _getFacePx(self) - P000 = V*Px('fXm') - P100 = V*Px('fXp') + fP = _getFacePx(self) + P000 = V*fP('fXm') + P100 = V*fP('fXp') elif d == 2: - Pxx = _getFacePxx(self) - P000 = V*Pxx('fXm', 'fYm') - P100 = V*Pxx('fXp', 'fYm') - P010 = V*Pxx('fXm', 'fYp') - P110 = V*Pxx('fXp', 'fYp') + fP = _getFacePxx(self) + P000 = V*fP('fXm', 'fYm') + P100 = V*fP('fXp', 'fYm') + P010 = V*fP('fXm', 'fYp') + P110 = V*fP('fXp', 'fYp') elif d == 3: - Pxxx = _getFacePxxx(self) - P000 = V*Pxxx('fXm', 'fYm', 'fZm') - P100 = V*Pxxx('fXp', 'fYm', 'fZm') - P010 = V*Pxxx('fXm', 'fYp', 'fZm') - P110 = V*Pxxx('fXp', 'fYp', 'fZm') - P001 = V*Pxxx('fXm', 'fYm', 'fZp') - P101 = V*Pxxx('fXp', 'fYm', 'fZp') - P011 = V*Pxxx('fXm', 'fYp', 'fZp') - P111 = V*Pxxx('fXp', 'fYp', 'fZp') + fP = _getFacePxxx(self) + P000 = V*fP('fXm', 'fYm', 'fZm') + P100 = V*fP('fXp', 'fYm', 'fZm') + P010 = V*fP('fXm', 'fYp', 'fZm') + P110 = V*fP('fXp', 'fYp', 'fZm') + P001 = V*fP('fXm', 'fYm', 'fZp') + P101 = V*fP('fXp', 'fYm', 'fZp') + P011 = V*fP('fXm', 'fYp', 'fZp') + P111 = V*fP('fXp', 'fYp', 'fZp') Mu = makePropertyTensor(self, materialProperty) A = P000.T*Mu*P000 + P100.T*Mu*P100 @@ -70,7 +70,6 @@ class InnerProducts(object): if d == 1: raise NotImplementedError('getEdgeInnerProduct not implemented for 1D') - # We will multiply by V on each side to keep symmetry elif d == 2: eP = _getEdgePxx(self) P000 = V*eP('eX0', 'eY0') diff --git a/docs/api_DiffOps.rst b/docs/api_DiffOps.rst new file mode 100644 index 00000000..6ddde039 --- /dev/null +++ b/docs/api_DiffOps.rst @@ -0,0 +1,8 @@ +.. _api_DiffOps: + +Differential Operators +====================== + +.. automodule:: SimPEG.Mesh.DiffOperators + :members: + :undoc-members: diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst index df5f31d8..2a3443d7 100644 --- a/docs/api_Mesh.rst +++ b/docs/api_Mesh.rst @@ -8,13 +8,13 @@ Tensor Mesh :show-inheritance: :members: :undoc-members: - :inherited-members: Cylindrical 1D Mesh =================== .. automodule:: SimPEG.Mesh.Cyl1DMesh + :show-inheritance: :members: :undoc-members: @@ -26,7 +26,6 @@ Logically Orthogonal Mesh :show-inheritance: :members: :undoc-members: - :inherited-members: Base Mesh @@ -37,9 +36,3 @@ Base Mesh :undoc-members: -Differential Operators -====================== - -.. automodule:: SimPEG.Mesh.DiffOperators - :members: - :undoc-members: diff --git a/docs/api_Solver.rst b/docs/api_Solver.rst new file mode 100644 index 00000000..be009f97 --- /dev/null +++ b/docs/api_Solver.rst @@ -0,0 +1,10 @@ +.. _api_Solver: + + +Solver +****** + +.. automodule:: SimPEG.Solver + :members: + :undoc-members: + diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst index 641e6d60..c8dfe47b 100644 --- a/docs/api_Utils.rst +++ b/docs/api_Utils.rst @@ -1,14 +1,6 @@ .. _api_Utils: -Solver -****** - -.. automodule:: SimPEG.Solver - :members: - :undoc-members: - - Utilities ********* diff --git a/docs/index.rst b/docs/index.rst index 6a282e36..e8da6b25 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -11,12 +11,10 @@ The vision is to create a package for finite volume simulation and parameter est applications to geophysical imaging and subsurface flow. To enable these goals, this package has the following features: - - * is modular with respect to discretization, physics, optimization, and regularization - * is built with the (large-scale) inverse problem in mind - * provides a framework for geophysical and hydrogeologic problems - * supports 1D, 2D and 3D problems - +- is modular with respect to discretization, physics, optimization, and regularization +- is built with the (large-scale) inverse problem in mind +- provides a framework for geophysical and hydrogeologic problems +- supports 1D, 2D and 3D problems .. image:: simpeg-framework.png :width: 400 px @@ -30,6 +28,7 @@ Meshing & Operators :maxdepth: 2 api_Mesh + api_DiffOps api_InnerProducts Forward Problems @@ -49,6 +48,16 @@ Inversion api_Inverse api_Parameters +Utility Codes +************* + +.. toctree:: + :maxdepth: 2 + + api_Solver + api_Utils + + Testing SimPEG ************** @@ -63,21 +72,6 @@ Testing SimPEG :alt: Master Branch :align: center -* Develop Branch - .. image:: https://travis-ci.org/simpeg/simpeg.png?branch*develop - :target: https://travis-ci.org/simpeg/simpeg - :alt: Develop Branch - :align: center - - -Utility Codes -************* - -.. toctree:: - :maxdepth: 2 - - api_Utils - Project Index & Search ********************** From d0a734b17c142fa2d3dc8f6e6b7b6f2e844d27cd Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 26 Feb 2014 16:29:33 -0800 Subject: [PATCH 04/24] Added option to invert the material property. --- SimPEG/Mesh/InnerProducts.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index a797a4a9..87ec7094 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -1,5 +1,5 @@ from scipy import sparse as sp -from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor +from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor, invPropertyTensor import numpy as np @@ -10,13 +10,19 @@ class InnerProducts(object): 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, materialProperty=None, returnP=False): + def getFaceInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False): """ :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 :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ + if invertProperty: + materialProperty = invPropertyTensor(self, materialProperty) + + Mu = makePropertyTensor(self, materialProperty) + d = self.dim # We will multiply by sqrt on each side to keep symmetry V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol))) @@ -42,7 +48,6 @@ class InnerProducts(object): P011 = V*fP('fXm', 'fYp', 'fZp') P111 = V*fP('fXp', 'fYp', 'fZp') - Mu = makePropertyTensor(self, materialProperty) A = P000.T*Mu*P000 + P100.T*Mu*P100 P = [P000, P100] @@ -57,13 +62,19 @@ class InnerProducts(object): else: return A - def getEdgeInnerProduct(self, materialProperty=None, returnP=False): + def getEdgeInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False): """ :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 :rtype: scipy.csr_matrix :return: M, the inner product matrix (nE, nE) """ + if invertProperty: + materialProperty = invPropertyTensor(self, materialProperty) + + Mu = makePropertyTensor(self, materialProperty) + d = self.dim # We will multiply by sqrt on each side to keep symmetry V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol))) From 69b920ec1cb347ce55091a67d0a6553c7af96805 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 27 Feb 2014 13:31:45 -0800 Subject: [PATCH 05/24] setup.py --- .gitignore | 37 +++- ez_setup.py | 485 ++++++++++++++++++++++++++++++++++++++++++++++++++++ setup.py | 46 +++++ 3 files changed, 567 insertions(+), 1 deletion(-) create mode 100644 ez_setup.py create mode 100644 setup.py diff --git a/.gitignore b/.gitignore index 7678b15c..14e3cc80 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,40 @@ -*.pyc +*.py[cod] + +# C extensions *.so + +# Packages +*.egg +*.egg-info +dist +build +eggs +parts +bin +var +sdist +develop-eggs +.installed.cfg +lib +lib64 +__pycache__ + +# Installer logs +pip-log.txt + +# Unit test / coverage reports +.coverage +.tox +nosetests.xml + +# Translations +*.mo + +# Mr Developer +.mr.developer.cfg +.project +.pydevproject + *.sublime-project *.sublime-workspace docs/_build/ diff --git a/ez_setup.py b/ez_setup.py new file mode 100644 index 00000000..3ea2e667 --- /dev/null +++ b/ez_setup.py @@ -0,0 +1,485 @@ +#!python +"""Bootstrap distribute installation + +If you want to use setuptools in your package's setup.py, just include this +file in the same directory with it, and add this to the top of your setup.py:: + + from distribute_setup import use_setuptools + use_setuptools() + +If you want to require a specific version of setuptools, set a download +mirror, or use an alternate download directory, you can do so by supplying +the appropriate options to ``use_setuptools()``. + +This file can also be run as a script to install or upgrade setuptools. +""" +import os +import sys +import time +import fnmatch +import tempfile +import tarfile +from distutils import log + +try: + from site import USER_SITE +except ImportError: + USER_SITE = None + +try: + import subprocess + + def _python_cmd(*args): + args = (sys.executable,) + args + return subprocess.call(args) == 0 + +except ImportError: + # will be used for python 2.3 + def _python_cmd(*args): + args = (sys.executable,) + args + # quoting arguments if windows + if sys.platform == 'win32': + def quote(arg): + if ' ' in arg: + return '"%s"' % arg + return arg + args = [quote(arg) for arg in args] + return os.spawnl(os.P_WAIT, sys.executable, *args) == 0 + +DEFAULT_VERSION = "0.6.14" +DEFAULT_URL = "http://pypi.python.org/packages/source/d/distribute/" +SETUPTOOLS_FAKED_VERSION = "0.6c11" + +SETUPTOOLS_PKG_INFO = """\ +Metadata-Version: 1.0 +Name: setuptools +Version: %s +Summary: xxxx +Home-page: xxx +Author: xxx +Author-email: xxx +License: xxx +Description: xxx +""" % SETUPTOOLS_FAKED_VERSION + + +def _install(tarball): + # extracting the tarball + tmpdir = tempfile.mkdtemp() + log.warn('Extracting in %s', tmpdir) + old_wd = os.getcwd() + try: + os.chdir(tmpdir) + tar = tarfile.open(tarball) + _extractall(tar) + tar.close() + + # going in the directory + subdir = os.path.join(tmpdir, os.listdir(tmpdir)[0]) + os.chdir(subdir) + log.warn('Now working in %s', subdir) + + # installing + log.warn('Installing Distribute') + if not _python_cmd('setup.py', 'install'): + log.warn('Something went wrong during the installation.') + log.warn('See the error message above.') + finally: + os.chdir(old_wd) + + +def _build_egg(egg, tarball, to_dir): + # extracting the tarball + tmpdir = tempfile.mkdtemp() + log.warn('Extracting in %s', tmpdir) + old_wd = os.getcwd() + try: + os.chdir(tmpdir) + tar = tarfile.open(tarball) + _extractall(tar) + tar.close() + + # going in the directory + subdir = os.path.join(tmpdir, os.listdir(tmpdir)[0]) + os.chdir(subdir) + log.warn('Now working in %s', subdir) + + # building an egg + log.warn('Building a Distribute egg in %s', to_dir) + _python_cmd('setup.py', '-q', 'bdist_egg', '--dist-dir', to_dir) + + finally: + os.chdir(old_wd) + # returning the result + log.warn(egg) + if not os.path.exists(egg): + raise IOError('Could not build the egg.') + + +def _do_download(version, download_base, to_dir, download_delay): + egg = os.path.join(to_dir, 'distribute-%s-py%d.%d.egg' + % (version, sys.version_info[0], sys.version_info[1])) + if not os.path.exists(egg): + tarball = download_setuptools(version, download_base, + to_dir, download_delay) + _build_egg(egg, tarball, to_dir) + sys.path.insert(0, egg) + import setuptools + setuptools.bootstrap_install_from = egg + + +def use_setuptools(version=DEFAULT_VERSION, download_base=DEFAULT_URL, + to_dir=os.curdir, download_delay=15, no_fake=True): + # making sure we use the absolute path + to_dir = os.path.abspath(to_dir) + was_imported = 'pkg_resources' in sys.modules or \ + 'setuptools' in sys.modules + try: + try: + import pkg_resources + if not hasattr(pkg_resources, '_distribute'): + if not no_fake: + _fake_setuptools() + raise ImportError + except ImportError: + return _do_download(version, download_base, to_dir, download_delay) + try: + pkg_resources.require("distribute>="+version) + return + except pkg_resources.VersionConflict: + e = sys.exc_info()[1] + if was_imported: + sys.stderr.write( + "The required version of distribute (>=%s) is not available,\n" + "and can't be installed while this script is running. Please\n" + "install a more recent version first, using\n" + "'easy_install -U distribute'." + "\n\n(Currently using %r)\n" % (version, e.args[0])) + sys.exit(2) + else: + del pkg_resources, sys.modules['pkg_resources'] # reload ok + return _do_download(version, download_base, to_dir, + download_delay) + except pkg_resources.DistributionNotFound: + return _do_download(version, download_base, to_dir, + download_delay) + finally: + if not no_fake: + _create_fake_setuptools_pkg_info(to_dir) + +def download_setuptools(version=DEFAULT_VERSION, download_base=DEFAULT_URL, + to_dir=os.curdir, delay=15): + """Download distribute from a specified location and return its filename + + `version` should be a valid distribute version number that is available + as an egg for download under the `download_base` URL (which should end + with a '/'). `to_dir` is the directory where the egg will be downloaded. + `delay` is the number of seconds to pause before an actual download + attempt. + """ + # making sure we use the absolute path + to_dir = os.path.abspath(to_dir) + try: + from urllib.request import urlopen + except ImportError: + from urllib2 import urlopen + tgz_name = "distribute-%s.tar.gz" % version + url = download_base + tgz_name + saveto = os.path.join(to_dir, tgz_name) + src = dst = None + if not os.path.exists(saveto): # Avoid repeated downloads + try: + log.warn("Downloading %s", url) + src = urlopen(url) + # Read/write all in one block, so we don't create a corrupt file + # if the download is interrupted. + data = src.read() + dst = open(saveto, "wb") + dst.write(data) + finally: + if src: + src.close() + if dst: + dst.close() + return os.path.realpath(saveto) + +def _no_sandbox(function): + def __no_sandbox(*args, **kw): + try: + from setuptools.sandbox import DirectorySandbox + if not hasattr(DirectorySandbox, '_old'): + def violation(*args): + pass + DirectorySandbox._old = DirectorySandbox._violation + DirectorySandbox._violation = violation + patched = True + else: + patched = False + except ImportError: + patched = False + + try: + return function(*args, **kw) + finally: + if patched: + DirectorySandbox._violation = DirectorySandbox._old + del DirectorySandbox._old + + return __no_sandbox + +def _patch_file(path, content): + """Will backup the file then patch it""" + existing_content = open(path).read() + if existing_content == content: + # already patched + log.warn('Already patched.') + return False + log.warn('Patching...') + _rename_path(path) + f = open(path, 'w') + try: + f.write(content) + finally: + f.close() + return True + +_patch_file = _no_sandbox(_patch_file) + +def _same_content(path, content): + return open(path).read() == content + +def _rename_path(path): + new_name = path + '.OLD.%s' % time.time() + log.warn('Renaming %s into %s', path, new_name) + os.rename(path, new_name) + return new_name + +def _remove_flat_installation(placeholder): + if not os.path.isdir(placeholder): + log.warn('Unkown installation at %s', placeholder) + return False + found = False + for file in os.listdir(placeholder): + if fnmatch.fnmatch(file, 'setuptools*.egg-info'): + found = True + break + if not found: + log.warn('Could not locate setuptools*.egg-info') + return + + log.warn('Removing elements out of the way...') + pkg_info = os.path.join(placeholder, file) + if os.path.isdir(pkg_info): + patched = _patch_egg_dir(pkg_info) + else: + patched = _patch_file(pkg_info, SETUPTOOLS_PKG_INFO) + + if not patched: + log.warn('%s already patched.', pkg_info) + return False + # now let's move the files out of the way + for element in ('setuptools', 'pkg_resources.py', 'site.py'): + element = os.path.join(placeholder, element) + if os.path.exists(element): + _rename_path(element) + else: + log.warn('Could not find the %s element of the ' + 'Setuptools distribution', element) + return True + +_remove_flat_installation = _no_sandbox(_remove_flat_installation) + +def _after_install(dist): + log.warn('After install bootstrap.') + placeholder = dist.get_command_obj('install').install_purelib + _create_fake_setuptools_pkg_info(placeholder) + +def _create_fake_setuptools_pkg_info(placeholder): + if not placeholder or not os.path.exists(placeholder): + log.warn('Could not find the install location') + return + pyver = '%s.%s' % (sys.version_info[0], sys.version_info[1]) + setuptools_file = 'setuptools-%s-py%s.egg-info' % \ + (SETUPTOOLS_FAKED_VERSION, pyver) + pkg_info = os.path.join(placeholder, setuptools_file) + if os.path.exists(pkg_info): + log.warn('%s already exists', pkg_info) + return + + log.warn('Creating %s', pkg_info) + f = open(pkg_info, 'w') + try: + f.write(SETUPTOOLS_PKG_INFO) + finally: + f.close() + + pth_file = os.path.join(placeholder, 'setuptools.pth') + log.warn('Creating %s', pth_file) + f = open(pth_file, 'w') + try: + f.write(os.path.join(os.curdir, setuptools_file)) + finally: + f.close() + +_create_fake_setuptools_pkg_info = _no_sandbox(_create_fake_setuptools_pkg_info) + +def _patch_egg_dir(path): + # let's check if it's already patched + pkg_info = os.path.join(path, 'EGG-INFO', 'PKG-INFO') + if os.path.exists(pkg_info): + if _same_content(pkg_info, SETUPTOOLS_PKG_INFO): + log.warn('%s already patched.', pkg_info) + return False + _rename_path(path) + os.mkdir(path) + os.mkdir(os.path.join(path, 'EGG-INFO')) + pkg_info = os.path.join(path, 'EGG-INFO', 'PKG-INFO') + f = open(pkg_info, 'w') + try: + f.write(SETUPTOOLS_PKG_INFO) + finally: + f.close() + return True + +_patch_egg_dir = _no_sandbox(_patch_egg_dir) + +def _before_install(): + log.warn('Before install bootstrap.') + _fake_setuptools() + + +def _under_prefix(location): + if 'install' not in sys.argv: + return True + args = sys.argv[sys.argv.index('install')+1:] + for index, arg in enumerate(args): + for option in ('--root', '--prefix'): + if arg.startswith('%s=' % option): + top_dir = arg.split('root=')[-1] + return location.startswith(top_dir) + elif arg == option: + if len(args) > index: + top_dir = args[index+1] + return location.startswith(top_dir) + if arg == '--user' and USER_SITE is not None: + return location.startswith(USER_SITE) + return True + + +def _fake_setuptools(): + log.warn('Scanning installed packages') + try: + import pkg_resources + except ImportError: + # we're cool + log.warn('Setuptools or Distribute does not seem to be installed.') + return + ws = pkg_resources.working_set + try: + setuptools_dist = ws.find(pkg_resources.Requirement.parse('setuptools', + replacement=False)) + except TypeError: + # old distribute API + setuptools_dist = ws.find(pkg_resources.Requirement.parse('setuptools')) + + if setuptools_dist is None: + log.warn('No setuptools distribution found') + return + # detecting if it was already faked + setuptools_location = setuptools_dist.location + log.warn('Setuptools installation detected at %s', setuptools_location) + + # if --root or --preix was provided, and if + # setuptools is not located in them, we don't patch it + if not _under_prefix(setuptools_location): + log.warn('Not patching, --root or --prefix is installing Distribute' + ' in another location') + return + + # let's see if its an egg + if not setuptools_location.endswith('.egg'): + log.warn('Non-egg installation') + res = _remove_flat_installation(setuptools_location) + if not res: + return + else: + log.warn('Egg installation') + pkg_info = os.path.join(setuptools_location, 'EGG-INFO', 'PKG-INFO') + if (os.path.exists(pkg_info) and + _same_content(pkg_info, SETUPTOOLS_PKG_INFO)): + log.warn('Already patched.') + return + log.warn('Patching...') + # let's create a fake egg replacing setuptools one + res = _patch_egg_dir(setuptools_location) + if not res: + return + log.warn('Patched done.') + _relaunch() + + +def _relaunch(): + log.warn('Relaunching...') + # we have to relaunch the process + # pip marker to avoid a relaunch bug + if sys.argv[:3] == ['-c', 'install', '--single-version-externally-managed']: + sys.argv[0] = 'setup.py' + args = [sys.executable] + sys.argv + sys.exit(subprocess.call(args)) + + +def _extractall(self, path=".", members=None): + """Extract all members from the archive to the current working + directory and set owner, modification time and permissions on + directories afterwards. `path' specifies a different directory + to extract to. `members' is optional and must be a subset of the + list returned by getmembers(). + """ + import copy + import operator + from tarfile import ExtractError + directories = [] + + if members is None: + members = self + + for tarinfo in members: + if tarinfo.isdir(): + # Extract directories with a safe mode. + directories.append(tarinfo) + tarinfo = copy.copy(tarinfo) + tarinfo.mode = 448 # decimal for oct 0700 + self.extract(tarinfo, path) + + # Reverse sort directories. + if sys.version_info < (2, 4): + def sorter(dir1, dir2): + return cmp(dir1.name, dir2.name) + directories.sort(sorter) + directories.reverse() + else: + directories.sort(key=operator.attrgetter('name'), reverse=True) + + # Set correct owner, mtime and filemode on directories. + for tarinfo in directories: + dirpath = os.path.join(path, tarinfo.name) + try: + self.chown(tarinfo, dirpath) + self.utime(tarinfo, dirpath) + self.chmod(tarinfo, dirpath) + except ExtractError: + e = sys.exc_info()[1] + if self.errorlevel > 1: + raise + else: + self._dbg(1, "tarfile: %s" % e) + + +def main(argv, version=DEFAULT_VERSION): + """Install or upgrade setuptools and EasyInstall""" + tarball = download_setuptools() + _install(tarball) + + +if __name__ == '__main__': + main(sys.argv[1:]) diff --git a/setup.py b/setup.py new file mode 100644 index 00000000..ed65ab33 --- /dev/null +++ b/setup.py @@ -0,0 +1,46 @@ +#!/usr/bin/env python +"""SimPEG: Simulation and Parameter Estimation for Geophysics + +SimPEG is a python package for simulation and gradient based +parameter estimation in the context of geophysical applications. +""" + +import ez_setup +ez_setup.use_setuptools() +from setuptools import setup, find_packages + +CLASSIFIERS = [ +'Development Status :: 0.0.1 - Alpha', +'Intended Audience :: Science/Research', +'Intended Audience :: Developers', +'License :: MIT License', +'Programming Language :: Python', +'Topic :: Scientific/Engineering', +'Topic :: Scientific/Engineering :: Mathematics', +'Operating System :: Microsoft :: Windows', +'Operating System :: POSIX', +'Operating System :: Unix', +'Operating System :: MacOS', +] + +import os, os.path + +setup( + name = "SimPEG", + version = "0.1dev", + packages = find_packages(), + install_requires = ['numpy>=1.7', + 'scipy>=0.12', + 'matplotlib>=1.3', + ], + author = "Rowan Cockett", + author_email = "rowanc1@gmail.com", + description = "SimPEG: Simulation and Parameter Estimation for Geophysics", + license = "MIT", + keywords = "geophysics inverse problem", + url = "http://simeg.rtfd.org/", + download_url = "http://github.com/simpeg", + classifiers=CLASSIFIERS, + platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"], + use_2to3 = False +) From e08467af3ec16e382d8b68fd30c114c746323352 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 28 Feb 2014 15:19:00 -0800 Subject: [PATCH 06/24] tolerance on test --- SimPEG/Tests/test_TreeMesh.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/SimPEG/Tests/test_TreeMesh.py b/SimPEG/Tests/test_TreeMesh.py index f8120859..7f2bce78 100644 --- a/SimPEG/Tests/test_TreeMesh.py +++ b/SimPEG/Tests/test_TreeMesh.py @@ -4,6 +4,8 @@ import numpy as np import unittest import matplotlib.pyplot as plt +TOL = 1e-10 + class TestOcTreeObjects(unittest.TestCase): def setUp(self): @@ -493,10 +495,10 @@ class SimpleOctreeOperatorTests(unittest.TestCase): # self.assertTrue((self.tM2.edgeCurl - self.oM2.edgeCurl).toarray().sum() == 0) def test_InnerProducts(self): - self.assertTrue((self.tM.getFaceInnerProduct() - self.oM.getFaceInnerProduct()).toarray().sum() == 0) - self.assertTrue((self.tM2.getFaceInnerProduct() - self.oM2.getFaceInnerProduct()).toarray().sum() == 0) - self.assertTrue((self.tM2.getEdgeInnerProduct() - self.oM2.getEdgeInnerProduct()).toarray().sum() == 0) - self.assertTrue((self.tM.getEdgeInnerProduct() - self.oM.getEdgeInnerProduct()).toarray().sum() == 0) + self.assertTrue((self.tM.getFaceInnerProduct() - self.oM.getFaceInnerProduct()).toarray().sum() < TOL) + self.assertTrue((self.tM2.getFaceInnerProduct() - self.oM2.getFaceInnerProduct()).toarray().sum() < TOL) + self.assertTrue((self.tM2.getEdgeInnerProduct() - self.oM2.getEdgeInnerProduct()).toarray().sum() < TOL) + self.assertTrue((self.tM.getEdgeInnerProduct() - self.oM.getEdgeInnerProduct()).toarray().sum() < TOL) if __name__ == '__main__': From 1e5053085d269e2c9bfe8a10625627943b76f02d Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 28 Feb 2014 15:19:57 -0800 Subject: [PATCH 07/24] changed matutils to handle vector material properties not just columns --- SimPEG/Tests/test_utils.py | 1 + SimPEG/Utils/matutils.py | 31 ++++++++++++++++++++----------- 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index 429c797c..e69c8473 100644 --- a/SimPEG/Tests/test_utils.py +++ b/SimPEG/Tests/test_utils.py @@ -71,6 +71,7 @@ class TestSequenceFunctions(unittest.TestCase): self.assertTrue(np.all(sub2ind(x.shape, [4,0]) == [4])) self.assertTrue(np.all(sub2ind(x.shape, [0,1]) == [5])) self.assertTrue(np.all(sub2ind(x.shape, [4,1]) == [9])) + self.assertTrue(np.all(sub2ind(x.shape, [[4,1]]) == [9])) self.assertTrue(np.all(sub2ind(x.shape, [[0,0],[4,0],[0,1],[4,1]]) == [0,4,5,9])) def test_ind2sub(self): diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index 9e350957..c15f6234 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -146,7 +146,7 @@ def sub2ind(shape, subs): """From the given shape, returns the index of the given subscript""" if type(subs) is not np.ndarray: subs = np.array(subs) - if subs.size == len(shape): + if len(subs.shape) == 1: subs = subs[np.newaxis,:] assert subs.shape[1] == len(shape), 'Indexing must be done as a column vectors. e.g. [[3,6],[6,2],...]' inds = np.ravel_multi_index(subs.T, shape, order='F') @@ -261,9 +261,11 @@ def makePropertyTensor(M, sigma): if sigma.size == M.nC: # Isotropic! sigma = mkvc(sigma) # ensure it is a vector. Sigma = sdiag(np.r_[sigma, sigma]) - elif sigma.shape[1] == 2: # Diagonal tensor + elif sigma.size == M.nC*2: # Diagonal tensor + sigma = sigma.reshape((M.nC,2), order='F') Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]]) - elif sigma.shape[1] == 3: # Fully anisotropic + elif sigma.size == M.nC*3: # Fully anisotropic + sigma = sigma.reshape((M.nC,3), order='F') row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2]))) row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1]))) Sigma = sp.vstack((row1, row2)) @@ -273,15 +275,18 @@ def makePropertyTensor(M, sigma): if sigma.size == M.nC: # Isotropic! sigma = mkvc(sigma) # ensure it is a vector. Sigma = sdiag(np.r_[sigma, sigma, sigma]) - elif sigma.shape[1] == 3: # Diagonal tensor + elif sigma.size == M.nC*3: # Diagonal tensor + sigma = sigma.reshape((M.nC,3), order='F') Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]]) - elif sigma.shape[1] == 6: # Fully anisotropic + elif sigma.size == M.nC*6: # Fully anisotropic + sigma = sigma.reshape((M.nC,6), order='F') row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4]))) row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5]))) row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2]))) Sigma = sp.vstack((row1, row2, row3)) else: raise Exception('Unexpected shape of sigma') + return Sigma @@ -296,27 +301,31 @@ def invPropertyTensor(M, tensor, returnMatrix=False): T = 1./mkvc(tensor) # ensure it is a vector. elif M.dim == 2: - if tensor.shape[1] == 2: # Diagonal tensor + if tensor.size == M.nC*2: # Diagonal tensor T = 1./tensor - elif tensor.shape[1] == 3: # Fully anisotropic + elif tensor.size == M.nC*3: # Fully anisotropic + tensor = tensor.reshape((M.nC,3), order='F') B = inv2X2BlockDiagonal(tensor[:,0], tensor[:,2], tensor[:,2], tensor[:,1], returnMatrix=False) b11, b12, b21, b22 = B - T = np.c_[b11, b22, b12] + T = np.r_[b11, b22, b12] elif M.dim == 3: - if tensor.shape[1] == 3: # Diagonal tensor + if tensor.size == M.nC*3: # Diagonal tensor T = 1./tensor - elif tensor.shape[1] == 6: # Fully anisotropic + elif tensor.size == M.nC*6: # Fully anisotropic + tensor = tensor.reshape((M.nC,6), order='F') B = inv3X3BlockDiagonal(tensor[:,0], tensor[:,3], tensor[:,4], tensor[:,3], tensor[:,1], tensor[:,5], tensor[:,4], tensor[:,5], tensor[:,2], returnMatrix=False) b11, b12, b13, b21, b22, b23, b31, b32, b33 = B - T = np.c_[b11, b22, b33, b12, b13, b23] + T = np.r_[b11, b22, b33, b12, b13, b23] if T is None: raise Exception('Unexpected shape of tensor') + if returnMatrix: return makePropertyTensor(M, T) + return T From b5c7b112650bc88247feb3708ae6e1f9e49d3d75 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 28 Feb 2014 15:21:39 -0800 Subject: [PATCH 08/24] updates to face innerProducts Increase Speed Add derivatives Add tests --- SimPEG/Mesh/InnerProducts.py | 39 +++++++++++++++++++++- SimPEG/Mesh/TensorMesh.py | 44 +++++++++++++++++++++++++ SimPEG/Tests/test_massMatrices.py | 2 +- SimPEG/Tests/test_massMatricesDerivs.py | 40 ++++++++++++++++++++++ docs/api_InnerProducts.rst | 37 +++++++++++++++++++++ 5 files changed, 160 insertions(+), 2 deletions(-) create mode 100644 SimPEG/Tests/test_massMatricesDerivs.py diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 87ec7094..034cdd01 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -10,7 +10,8 @@ class InnerProducts(object): 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, materialProperty=None, returnP=False, invertProperty=False): + def getFaceInnerProduct(self, materialProperty=None, returnP=False, + invertProperty=False): """ :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 @@ -18,6 +19,14 @@ class InnerProducts(object): :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ + fast = None + + if returnP is False and hasattr(self, '_fastFaceInnerProduct'): + fast = self._fastFaceInnerProduct(materialProperty=materialProperty, invertProperty=invertProperty) + + if fast is not None: + return fast + if invertProperty: materialProperty = invPropertyTensor(self, materialProperty) @@ -62,6 +71,34 @@ class InnerProducts(object): else: return A + def getFaceInnerProductDeriv(self, materialProperty=None, P=None): + """ + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + + fast = None + + if hasattr(self, '_fastFaceInnerProductDeriv'): + fast = self._fastFaceInnerProductDeriv(materialProperty=materialProperty) + + if fast is not None: + return fast + + raise NotImplementedError('Derivatives for the material property specified are not yet implemented.') + + if P is None: + M, P = getFaceInnerProduct(self, materialProperty=materialProperty, returnP=True) + + d = self.dim + + if d == 1: + P[0].T * sp.identity(n) * P[0] + + + + def getEdgeInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 8a80ffea..f55c6df6 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -39,6 +39,7 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): def __init__(self, h_in, x0=None): assert type(h_in) is list, 'h_in must be a list' + assert len(h_in) in [1,2,3], 'h_in must be of dimension 1, 2, or 3' h = range(len(h_in)) for i, h_i in enumerate(h_in): if type(h_i) in [int, long, float]: @@ -491,6 +492,49 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): indzu = (self.gridCC[:,2]==max(self.gridCC[:,2])) return indxd, indxu, indyd, indyu, indzd, indzu + def _fastFaceInnerProduct(self, materialProperty=None, invertProperty=False): + """ + Fast version of getFaceInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :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 + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + if materialProperty is None: + materialProperty = np.ones(self.nC) + + if materialProperty.size == self.nC: + if invertProperty: + materialProperty = 1./materialProperty + Av = self.aveF2CC + V = Utils.sdiag(self.vol) + return self.dim * Utils.sdiag(Av.T * V * materialProperty) + if materialProperty.size == self.nC*self.dim: + if invertProperty: + materialProperty = 1./materialProperty + Av = self.aveF2CCV + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty)) + + + + def _fastFaceInnerProductDeriv(self, materialProperty=None): + """ + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + if materialProperty is None or materialProperty.size == self.nC: + Av = self.aveF2CC + return self.dim * Av.T * Utils.sdiag(self.vol) + if materialProperty.size == self.nC*self.dim: # anisotropic + Av = self.aveF2CCV + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + return Av.T * V + if __name__ == '__main__': print('Welcome to tensor mesh!') diff --git a/SimPEG/Tests/test_massMatrices.py b/SimPEG/Tests/test_massMatrices.py index 0db8c007..200d2884 100644 --- a/SimPEG/Tests/test_massMatrices.py +++ b/SimPEG/Tests/test_massMatrices.py @@ -30,7 +30,7 @@ class TestInnerProducts(OrderTest): sigma = np.c_[call(sigma1, Gc)] analytic = 647./360 # Found using sympy. elif self.sigmaTest == 3: - sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] + sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] analytic = 37./12 # Found using sympy. elif self.sigmaTest == 6: sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc), diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py new file mode 100644 index 00000000..f5484568 --- /dev/null +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -0,0 +1,40 @@ +import numpy as np +import unittest +from SimPEG import * +from TestUtils import checkDerivative + + +class TestInnerProductsDerivs(unittest.TestCase): + def setUp(self): + pass + + def test_FaceIP_derivs_isotropic(self): + for d in range(3): + mesh = Mesh.TensorMesh([10,5,4][d:]) + M,Ps = mesh.getFaceInnerProduct(returnP=True) + v = np.random.rand(mesh.nF) + def fun(sig): + M = mesh.getFaceInnerProduct(sig) + Md = mesh.getFaceInnerProductDeriv(sig) + return M*v, Utils.sdiag(v)*Md + sig = np.random.rand(mesh.nC) + passed = checkDerivative(fun, sig, plotIt=False) + self.assertTrue(passed) + + def test_FaceIP_derivs_anisotropic(self): + for d in range(3): + mesh = Mesh.TensorMesh([10,5,4][d:]) + M,Ps = mesh.getFaceInnerProduct(returnP=True) + v = np.random.rand(mesh.nF) + def fun(sig): + M = mesh.getFaceInnerProduct(sig) + Md = mesh.getFaceInnerProductDeriv(sig) + return M*v, Utils.sdiag(v)*Md + sig = np.random.rand(mesh.nC*mesh.dim) + passed = checkDerivative(fun, sig, plotIt=False) + self.assertTrue(passed) + + + +if __name__ == '__main__': + unittest.main() diff --git a/docs/api_InnerProducts.rst b/docs/api_InnerProducts.rst index dec9c4c8..06c0b2ed 100644 --- a/docs/api_InnerProducts.rst +++ b/docs/api_InnerProducts.rst @@ -131,6 +131,7 @@ If ``returnP=True`` is requested in any of these methods the projection matrices The derivation for ``edgeInnerProducts`` is exactly the same, however, when we approximate the integral using the fields around each node, the projection matrices look a bit different because we have 12 edges in 3D instead of just 6 faces. The interface to the code is exactly the same. + Defining Tensor Properties -------------------------- @@ -158,6 +159,42 @@ Depending on the number of columns (either 1, 3, or 6) of mu, the material prope \vec{\mu} = \left[\begin{matrix} \mu_{1} & \mu_{3} \\ \mu_{3} & \mu_{2} \end{matrix}\right] +Structure of Matrices +--------------------- + +Both the isotropic, and anisotropic material properties result in a diagonal mass matrix. +Which is nice and easy to invert if necessary, however, in the fully anisotropic case which is not aligned with the grid, the matrix is not diagonal. This can be seen for a 3D mesh in the figure below. + +.. plot:: + + from SimPEG import * + mesh = Mesh.TensorMesh([10,50,3]) + m1 = np.random.rand(mesh.nC) + m2 = np.random.rand(mesh.nC,3) + m3 = np.random.rand(mesh.nC,6) + M = range(3) + M[0] = mesh.getFaceInnerProduct(m1) + M[1] = mesh.getFaceInnerProduct(m2) + M[2] = mesh.getFaceInnerProduct(m3) + plt.figure(figsize=(13,5)) + for i, lab in enumerate(['Isotropic','Anisotropic','Tensor']): + plt.subplot(131 + i) + plt.spy(M[i],ms=0.5,color='k') + plt.tick_params(axis='both',which='both',labeltop='off',labelleft='off') + plt.title(lab + ' Material Property') + plt.show() + + +Taking Derivatives +------------------ + +TODO: Take the derivatives of the tensors. + +.. math:: + + \left[\begin{smallmatrix}0.5 \sigma_{1} & 0 & 0.25 \sigma_{3} & 0.25 \sigma_{3}\\0 & 0.5 \sigma_{1} & 0.25 \sigma_{3} & 0.25 \sigma_{3}\\0.25 \sigma_{3} & 0.25 \sigma_{3} & 0.5 \sigma_{2} & 0\\0.25 \sigma_{3} & 0.25 \sigma_{3} & 0 & 0.5 \sigma_{2}\end{smallmatrix}\right] + + The API ------- From 3eb8224e28609f0abc46d07a231e8895f17d2ffc Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 1 Mar 2014 06:54:30 -0800 Subject: [PATCH 09/24] Edge inner products are derivatives w/ testing. --- SimPEG/Mesh/InnerProducts.py | 30 +++++++++++-- SimPEG/Mesh/TensorMesh.py | 57 ++++++++++++++++++++++--- SimPEG/Tests/test_massMatricesDerivs.py | 15 ++++++- 3 files changed, 92 insertions(+), 10 deletions(-) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 034cdd01..266661a7 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -77,7 +77,6 @@ class InnerProducts(object): :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ - fast = None if hasattr(self, '_fastFaceInnerProductDeriv'): @@ -97,8 +96,6 @@ class InnerProducts(object): P[0].T * sp.identity(n) * P[0] - - def getEdgeInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) @@ -107,6 +104,14 @@ class InnerProducts(object): :rtype: scipy.csr_matrix :return: M, the inner product matrix (nE, nE) """ + fast = None + + if returnP is False and hasattr(self, '_fastEdgeInnerProduct'): + fast = self._fastEdgeInnerProduct(materialProperty=materialProperty, invertProperty=invertProperty) + + if fast is not None: + return fast + if invertProperty: materialProperty = invPropertyTensor(self, materialProperty) @@ -146,6 +151,25 @@ class InnerProducts(object): else: return A + + def getEdgeInnerProductDeriv(self, materialProperty=None, P=None): + """ + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + + fast = None + + if hasattr(self, '_fastEdgeInnerProductDeriv'): + fast = self._fastEdgeInnerProductDeriv(materialProperty=materialProperty) + + if fast is not None: + return fast + + raise NotImplementedError('Derivatives for the material property specified are not yet implemented.') + + # ------------------------ Geometries ------------------------------ # # diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index f55c6df6..5fd53344 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -503,35 +503,82 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ + return self._fastInnerProduct('F', materialProperty=materialProperty, invertProperty=invertProperty) + + + def _fastEdgeInnerProduct(self, materialProperty=None, invertProperty=False): + """ + Fast version of getEdgeInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :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 + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nE, nE) + """ + return self._fastInnerProduct('E', materialProperty=materialProperty, invertProperty=invertProperty) + + + def _fastInnerProduct(self, AvType, materialProperty=None, invertProperty=False): + """ + Fast version of getFaceInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param str AvType: 'E' or 'F' + :param bool returnP: returns the projection matrices + :param bool invertProperty: inverts the material property + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ if materialProperty is None: materialProperty = np.ones(self.nC) if materialProperty.size == self.nC: if invertProperty: materialProperty = 1./materialProperty - Av = self.aveF2CC + Av = getattr(self, 'ave'+AvType+'2CC') V = Utils.sdiag(self.vol) return self.dim * Utils.sdiag(Av.T * V * materialProperty) if materialProperty.size == self.nC*self.dim: if invertProperty: materialProperty = 1./materialProperty - Av = self.aveF2CCV + Av = getattr(self, 'ave'+AvType+'2CCV') V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty)) - def _fastFaceInnerProductDeriv(self, materialProperty=None): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ + return self._fastInnerProductDeriv('F', materialProperty=materialProperty) + + + def _fastEdgeInnerProductDeriv(self, materialProperty=None): + """ + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nE, nE) + """ + return self._fastInnerProductDeriv('E', materialProperty=materialProperty) + + + def _fastInnerProductDeriv(self, AvType, materialProperty=None): + """ + :param str AvType: 'E' or 'F' + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ if materialProperty is None or materialProperty.size == self.nC: - Av = self.aveF2CC + Av = getattr(self, 'ave'+AvType+'2CC') return self.dim * Av.T * Utils.sdiag(self.vol) if materialProperty.size == self.nC*self.dim: # anisotropic - Av = self.aveF2CCV + Av = getattr(self, 'ave'+AvType+'2CCV') V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) return Av.T * V diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py index f5484568..80f6bb83 100644 --- a/SimPEG/Tests/test_massMatricesDerivs.py +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -11,7 +11,6 @@ class TestInnerProductsDerivs(unittest.TestCase): def test_FaceIP_derivs_isotropic(self): for d in range(3): mesh = Mesh.TensorMesh([10,5,4][d:]) - M,Ps = mesh.getFaceInnerProduct(returnP=True) v = np.random.rand(mesh.nF) def fun(sig): M = mesh.getFaceInnerProduct(sig) @@ -21,10 +20,22 @@ class TestInnerProductsDerivs(unittest.TestCase): passed = checkDerivative(fun, sig, plotIt=False) self.assertTrue(passed) + + def test_EdgeIP_derivs_isotropic(self): + for h in [[10,5],[10,5,4]]: + mesh = Mesh.TensorMesh(h) + v = np.random.rand(mesh.nE) + def fun(sig): + M = mesh.getEdgeInnerProduct(sig) + Md = mesh.getEdgeInnerProductDeriv(sig) + return M*v, Utils.sdiag(v)*Md + sig = np.random.rand(mesh.nC) + passed = checkDerivative(fun, sig, plotIt=False) + self.assertTrue(passed) + def test_FaceIP_derivs_anisotropic(self): for d in range(3): mesh = Mesh.TensorMesh([10,5,4][d:]) - M,Ps = mesh.getFaceInnerProduct(returnP=True) v = np.random.rand(mesh.nF) def fun(sig): M = mesh.getFaceInnerProduct(sig) From eb91dd0026a697c38eedcc5be5fb530de9163364 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 2 Mar 2014 09:35:18 -0800 Subject: [PATCH 10/24] Updated linear inversion --- Tutorials/Linear.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Tutorials/Linear.py b/Tutorials/Linear.py index e18750c1..7d3aa805 100644 --- a/Tutorials/Linear.py +++ b/Tutorials/Linear.py @@ -12,10 +12,10 @@ class LinearProblem(Problem.BaseProblem): def fields(self, m, u=None): return self.G.dot(m) - def J(self, m, v, u=None): + def Jvec(self, m, v, u=None): return self.G.dot(v) - def Jt(self, m, v, u=None): + def Jtvec(self, m, v, u=None): return self.G.T.dot(v) From 27ae74d34aba0b314dc5a11a1805aa59ec42acca Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 2 Mar 2014 09:51:04 -0800 Subject: [PATCH 11/24] bug fix for windows compatability on solver. --- SimPEG/Solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index 6c47be27..68f36798 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -61,7 +61,7 @@ class Solver(object): warnings.warn("You should provide a preconditioner, M.", UserWarning) return M = options['M'] - if type(M) is sp.linalg.LinearOperator: + if isinstance(M, sp.linalg.LinearOperator): return PreconditionerList = ['J','GS'] if type(M) is str: From 908b9ab4cca9ce989f2dbb67bb2418e042e7cd5e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 2 Mar 2014 18:37:38 -0800 Subject: [PATCH 12/24] bug fix n --> vnC --- SimPEG/Mesh/DiffOperators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 03536f88..05c79b73 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -206,7 +206,7 @@ class DiffOperators(object): if(self.dim < 3): return None if(self._faceDivz is None): # The number of cell centers in each direction - n = self.n + n = self.vnC # Compute faceDivergence operator on faces D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0])) # Compute areas of cell faces & volumes From a01ffba4940f9078bc2818e81132bd7a3581dfcf Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 2 Mar 2014 18:39:31 -0800 Subject: [PATCH 13/24] make modelPair be OK in problem if first ComboModel is the right one. --- SimPEG/Problem.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 414d0c33..a509f6cd 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -44,7 +44,9 @@ class BaseProblem(object): def __init__(self, mesh, model, **kwargs): Utils.setKwargs(self, **kwargs) self.mesh = mesh - assert isinstance(model, self.modelPair), "Model object must be an instance of a %s class."%(self.modelPair.__name__) + assert (isinstance(model, self.modelPair) or + isinstance(model, Model.ComboModel) and isinstance(model.models[0], self.modelPair) + ), "Model object must be an instance of a %s class."%(self.modelPair.__name__) self.model = model @property From 802094096ec536a69e10b151a9e950d516cf7fc7 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 2 Mar 2014 18:46:56 -0800 Subject: [PATCH 14/24] change n --> vnC in cellGradz --- SimPEG/Mesh/DiffOperators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 05c79b73..a7a7b0cb 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -407,7 +407,7 @@ class DiffOperators(object): if self.dim < 3: return None if getattr(self, '_cellGradz', None) is None: BC = ['neumann', 'neumann'] - n = self.n + n = self.vnC G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0])) # Compute areas of cell faces & volumes V = self.aveCC2F*self.vol From c1d37ef8ffe3ada2c06cf166ce12365980e4c1b6 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 2 Mar 2014 20:15:29 -0800 Subject: [PATCH 15/24] change modelObj2Deriv from (self) --> (self, m, v=None) Had to change location of _bfgsH0 calculation, which now actually works. --- SimPEG/Inversion.py | 6 ------ SimPEG/ObjFunction.py | 2 +- SimPEG/Optimization.py | 14 ++++++++++---- SimPEG/Parameters.py | 2 +- SimPEG/Regularization.py | 17 ++++++++++++++--- 5 files changed, 26 insertions(+), 15 deletions(-) diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py index 64445ab2..297270ef 100644 --- a/SimPEG/Inversion.py +++ b/SimPEG/Inversion.py @@ -32,12 +32,6 @@ class BaseInversion(object): self.opt.printers.insert(2,IterationPrinters.phi_d) self.opt.printers.insert(3,IterationPrinters.phi_m) - if not hasattr(opt, '_bfgsH0') and hasattr(opt, 'bfgsH0'): # Check if it has been set by the user and the default is not being used. - #TODO: I don't think that this if statement is working... - print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.' - opt.bfgsH0 = SimPEG.Solver(objFunc.reg.modelObj2Deriv()) - - #TODO: Move this to the data class? @property def phi_d_target(self): diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py index 47354b4b..2b1912b8 100644 --- a/SimPEG/ObjFunction.py +++ b/SimPEG/ObjFunction.py @@ -97,7 +97,7 @@ class BaseObjFunction(object): if return_H: def H_fun(v): phi_d2Deriv = self.dataObj2Deriv(m, v, u=u) - phi_m2Deriv = self.reg.modelObj2Deriv()*v + phi_m2Deriv = self.reg.modelObj2Deriv(m, v=v) return phi_d2Deriv + self.beta * phi_m2Deriv diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 8e08250b..8bb5c95c 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -648,10 +648,16 @@ class BFGS(Minimize, Remember): Must be a SimPEG.Solver """ - _bfgsH0 = getattr(self,'_bfgsH0',None) - if _bfgsH0 is None: - return Solver(sp.identity(self.xc.size).tocsc(), flag='D') - return _bfgsH0 + if getattr(self,'_bfgsH0',None) is None: + # Check if it has been set by the user and the default is not being used. + if self.parent is None: + self._bfgsH0 = Solver(sp.identity(self.xc.size).tocsc(), flag='D') + else: + print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.' + objFunc = self.parent.objFunc + self._bfgsH0 = Solver(objFunc.reg.modelObj2Deriv(objFunc.m_current)) + return self._bfgsH0 + @bfgsH0.setter def bfgsH0(self, value): assert type(value) is Solver, 'bfgsH0 must be a SimPEG.Solver' diff --git a/SimPEG/Parameters.py b/SimPEG/Parameters.py index 9b12d8d1..5a6760e6 100644 --- a/SimPEG/Parameters.py +++ b/SimPEG/Parameters.py @@ -141,7 +141,7 @@ class BetaEstimate(Parameter): x0 = np.random.rand(*m.shape) t = x0.dot(objFunc.dataObj2Deriv(m,x0,u=u)) - b = x0.dot(objFunc.reg.modelObj2Deriv()*x0) + b = x0.dot(objFunc.reg.modelObj2Deriv(m, v=x0)) return self.beta0_ratio*(t/b) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 4fbf5d4b..8dc32b79 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -79,12 +79,18 @@ class BaseRegularization(object): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - return self.W.T * ( self.W * self.model.transform(m - self.mref) ) + mTd = self.model.transformDeriv(m - self.mref) + return mTd.T * ( self.W.T * ( self.W * self.model.transform(m - self.mref) ) ) @Utils.timeIt - def modelObj2Deriv(self): + def modelObj2Deriv(self, m, v=None): """ + :param numpy.array m: geophysical model + :param numpy.array v: vector to multiply + :rtype: scipy.sparse.csr_matrix or numpy.ndarray + :return: WtW or WtW*v + The regularization is: .. math:: @@ -98,7 +104,12 @@ class BaseRegularization(object): R(m) = \mathbf{W^\\top W} """ - return self.W.T * self.W + mTd = self.model.transformDeriv(m - self.mref) + if v is None: + return mTd.T * self.W.T * self.W * mTd + + return mTd.T * ( self.W.T * ( self.W * ( mTd * v) ) ) + From 8eaff94e48bb752eba62d2d501a3393b1946a16c Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 2 Mar 2014 23:34:31 -0800 Subject: [PATCH 16/24] Documentation for taking the derivatives of innerProducts --- docs/api_InnerProducts.rst | 72 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 70 insertions(+), 2 deletions(-) diff --git a/docs/api_InnerProducts.rst b/docs/api_InnerProducts.rst index 06c0b2ed..c3fbf4d6 100644 --- a/docs/api_InnerProducts.rst +++ b/docs/api_InnerProducts.rst @@ -188,12 +188,80 @@ Which is nice and easy to invert if necessary, however, in the fully anisotropic Taking Derivatives ------------------ -TODO: Take the derivatives of the tensors. +We will take the derivative of the fully anisotropic tensor for a 3D mesh, the other cases are easier and will not be discussed here. Let us start with one part of the sum which makes up \\\(\\mathbf{M}^f_\\Sigma\\\) and take the derivative when this is multiplied by some vector \\\(\\mathbf{v}\\\): .. math:: - \left[\begin{smallmatrix}0.5 \sigma_{1} & 0 & 0.25 \sigma_{3} & 0.25 \sigma_{3}\\0 & 0.5 \sigma_{1} & 0.25 \sigma_{3} & 0.25 \sigma_{3}\\0.25 \sigma_{3} & 0.25 \sigma_{3} & 0.5 \sigma_{2} & 0\\0.25 \sigma_{3} & 0.25 \sigma_{3} & 0 & 0.5 \sigma_{2}\end{smallmatrix}\right] + \mathbf{P}^\top \boldsymbol{\Sigma} \mathbf{Pv} +Here we will let \\\( \\mathbf{Pv} = \\mathbf{y} \\\) and \\\(\\mathbf{y}\\\) will have the form: + +.. math:: + + \mathbf{y} = \mathbf{Pv} = + \left[ + \begin{matrix} + \mathbf{y}_1\\ + \mathbf{y}_2\\ + \mathbf{y}_3\\ + \end{matrix} + \right] + +.. math:: + + \mathbf{P}^\top\Sigma\mathbf{y} = + \mathbf{P}^\top + \left[\begin{matrix} + \boldsymbol{\sigma}_{1} & \boldsymbol{\sigma}_{4} & \boldsymbol{\sigma}_{5} \\ + \boldsymbol{\sigma}_{4} & \boldsymbol{\sigma}_{2} & \boldsymbol{\sigma}_{6} \\ + \boldsymbol{\sigma}_{5} & \boldsymbol{\sigma}_{6} & \boldsymbol{\sigma}_{3} + \end{matrix}\right] + \left[ + \begin{matrix} + \mathbf{y}_1\\ + \mathbf{y}_2\\ + \mathbf{y}_3\\ + \end{matrix} + \right] + = + \mathbf{P}^\top + \left[ + \begin{matrix} + \boldsymbol{\sigma}_{1}\circ \mathbf{y}_1 + \boldsymbol{\sigma}_{4}\circ \mathbf{y}_2 + \boldsymbol{\sigma}_{5}\circ \mathbf{y}_3\\ + \boldsymbol{\sigma}_{4}\circ \mathbf{y}_1 + \boldsymbol{\sigma}_{2}\circ \mathbf{y}_2 + \boldsymbol{\sigma}_{6}\circ \mathbf{y}_3\\ + \boldsymbol{\sigma}_{5}\circ \mathbf{y}_1 + \boldsymbol{\sigma}_{6}\circ \mathbf{y}_2 + \boldsymbol{\sigma}_{3}\circ \mathbf{y}_3\\ + \end{matrix} + \right] + +Now it is easy to take the derivative with respect to any one of the parameters, for example, \\\(\\frac{\\partial}{\\partial\\boldsymbol{\\sigma}_1}\\\) + +.. math:: + \frac{\partial}{\partial \boldsymbol{\sigma}_1}\left(\mathbf{P}^\top\Sigma\mathbf{y}\right) + = + \mathbf{P}^\top + \left[ + \begin{matrix} + \text{diag}(\mathbf{y}_1)\\ + 0\\ + 0\\ + \end{matrix} + \right] + +Whereas \\\(\\frac{\\partial}{\\partial\\boldsymbol{\\sigma}_4}\\\), for example, is: + +.. math:: + \frac{\partial}{\partial \boldsymbol{\sigma}_4}\left(\mathbf{P}^\top\Sigma\mathbf{y}\right) + = + \mathbf{P}^\top + \left[ + \begin{matrix} + \text{diag}(\mathbf{y}_2)\\ + \text{diag}(\mathbf{y}_1)\\ + 0\\ + \end{matrix} + \right] + +These are computed for each of the 8 projections, horizontally concatenated, and returned. The API ------- From 3dd54bc0e3733301e69d54d89c9296fafcf98902 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 3 Mar 2014 11:01:12 -0800 Subject: [PATCH 17/24] Full anisotropic tensor derivatives. --- SimPEG/Mesh/InnerProducts.py | 123 ++++++++++++++++++++---- SimPEG/Mesh/TensorMesh.py | 19 ++-- SimPEG/Tests/test_massMatricesDerivs.py | 107 ++++++++++++++------- 3 files changed, 186 insertions(+), 63 deletions(-) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 266661a7..bb77d5e5 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -1,5 +1,5 @@ from scipy import sparse as sp -from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor, invPropertyTensor +from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor, invPropertyTensor, spzeros import numpy as np @@ -11,7 +11,7 @@ class InnerProducts(object): 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, materialProperty=None, returnP=False, - invertProperty=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 @@ -21,7 +21,7 @@ class InnerProducts(object): """ fast = None - if returnP is False and hasattr(self, '_fastFaceInnerProduct'): + if returnP is False and hasattr(self, '_fastFaceInnerProduct') and doFast: fast = self._fastFaceInnerProduct(materialProperty=materialProperty, invertProperty=invertProperty) if fast is not None: @@ -71,7 +71,7 @@ class InnerProducts(object): else: return A - def getFaceInnerProductDeriv(self, materialProperty=None, P=None): + def getFaceInnerProductDeriv(self, materialProperty=None, v=None, P=None, doFast=True): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :rtype: scipy.csr_matrix @@ -79,24 +79,18 @@ class InnerProducts(object): """ fast = None - if hasattr(self, '_fastFaceInnerProductDeriv'): - fast = self._fastFaceInnerProductDeriv(materialProperty=materialProperty) + if hasattr(self, '_fastFaceInnerProductDeriv') and doFast: + fast = self._fastFaceInnerProductDeriv(materialProperty=materialProperty, v=v) if fast is not None: return fast - raise NotImplementedError('Derivatives for the material property specified are not yet implemented.') - if P is None: - M, P = getFaceInnerProduct(self, materialProperty=materialProperty, returnP=True) + M, P = self.getFaceInnerProduct(materialProperty=materialProperty, returnP=True) - d = self.dim + return self._getInnerProductDeriv(materialProperty, v, P, self.nF) - if d == 1: - P[0].T * sp.identity(n) * P[0] - - - def getEdgeInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False): + 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 @@ -106,7 +100,7 @@ class InnerProducts(object): """ fast = None - if returnP is False and hasattr(self, '_fastEdgeInnerProduct'): + if returnP is False and hasattr(self, '_fastEdgeInnerProduct') and doFast: fast = self._fastEdgeInnerProduct(materialProperty=materialProperty, invertProperty=invertProperty) if fast is not None: @@ -151,8 +145,7 @@ class InnerProducts(object): else: return A - - def getEdgeInnerProductDeriv(self, materialProperty=None, P=None): + def getEdgeInnerProductDeriv(self, materialProperty=None, v=None, P=None, doFast=True): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :rtype: scipy.csr_matrix @@ -161,14 +154,102 @@ class InnerProducts(object): fast = None - if hasattr(self, '_fastEdgeInnerProductDeriv'): - fast = self._fastEdgeInnerProductDeriv(materialProperty=materialProperty) + if hasattr(self, '_fastEdgeInnerProductDeriv') and doFast: + fast = self._fastEdgeInnerProductDeriv(materialProperty=materialProperty, v=v) if fast is not None: return fast - raise NotImplementedError('Derivatives for the material property specified are not yet implemented.') + if P is None: + M, P = self.getEdgeInnerProduct(materialProperty=materialProperty, returnP=True) + return self._getInnerProductDeriv(materialProperty, v, P, self.nE) + + def _getInnerProductDeriv(self, materialProperty, v, P, n): + + if v is None: + raise Exception('v must be supplied for this implementation.') + + d = self.dim + Z = spzeros(self.nC, self.nC) + + if d == 1: + dMdm = spzeros(n, self.nC) + for i, p in enumerate(P): + dMdm = dMdm + p.T * sdiag( p * v ) + elif d == 2: + if materialProperty is None or materialProperty.size == self.nC: + dMdm = spzeros(n, self.nC) + for i, p in enumerate(P): + Y = p * v + y1 = Y[:self.nC] + y2 = Y[self.nC:] + dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 ))) + if materialProperty.size == self.nC*2: + dMdm1 = spzeros(n, self.nC) + dMdm2 = spzeros(n, self.nC) + for i, p in enumerate(P): + Y = p * v + y1 = Y[:self.nC] + y2 = Y[self.nC:] + dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z)) + dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ))) + dMdm = sp.hstack((dMdm1, dMdm2)) + if materialProperty.size == self.nC*3: + dMdm1 = spzeros(n, self.nC) + dMdm2 = spzeros(n, self.nC) + dMdm3 = spzeros(n, self.nC) + for i, p in enumerate(P): + Y = p * v + y1 = Y[:self.nC] + y2 = Y[self.nC:] + dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z)) + dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ))) + dMdm3 = dMdm3 + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 ))) + dMdm = sp.hstack((dMdm1, dMdm2, dMdm3)) + elif d == 3: + if materialProperty is None or materialProperty.size == self.nC: + dMdm = spzeros(n, self.nC) + for i, p in enumerate(P): + Y = p * v + y1 = Y[:self.nC] + y2 = Y[self.nC:self.nC*2] + y3 = Y[self.nC*2:] + dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 ), sdiag( y3 ))) + if materialProperty.size == self.nC*3: + dMdm1 = spzeros(n, self.nC) + dMdm2 = spzeros(n, self.nC) + dMdm3 = spzeros(n, self.nC) + for i, p in enumerate(P): + Y = p * v + y1 = Y[:self.nC] + y2 = Y[self.nC:self.nC*2] + y3 = Y[self.nC*2:] + dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z, Z)) + dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ), Z)) + dMdm3 = dMdm3 + p.T * sp.vstack(( Z, Z, sdiag( y3 ))) + dMdm = sp.hstack((dMdm1, dMdm2, dMdm3)) + if materialProperty.size == self.nC*6: + dMdm1 = spzeros(n, self.nC) + dMdm2 = spzeros(n, self.nC) + dMdm3 = spzeros(n, self.nC) + dMdm4 = spzeros(n, self.nC) + dMdm5 = spzeros(n, self.nC) + dMdm6 = spzeros(n, self.nC) + for i, p in enumerate(P): + Y = p * v + y1 = Y[:self.nC] + y2 = Y[self.nC:self.nC*2] + y3 = Y[self.nC*2:] + dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z, Z)) + dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ), Z)) + dMdm3 = dMdm3 + p.T * sp.vstack(( Z, Z, sdiag( y3 ))) + dMdm4 = dMdm4 + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 ), Z)) + dMdm5 = dMdm5 + p.T * sp.vstack(( sdiag( y3 ), Z, sdiag( y1 ))) + dMdm6 = dMdm6 + p.T * sp.vstack(( Z, sdiag( y3 ), sdiag( y2 ))) + dMdm = sp.hstack((dMdm1, dMdm2, dMdm3, dMdm4, dMdm5, dMdm6)) + + return dMdm # ------------------------ Geometries ------------------------------ # diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 5fd53344..e142d78d 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -549,25 +549,25 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty)) - def _fastFaceInnerProductDeriv(self, materialProperty=None): + def _fastFaceInnerProductDeriv(self, materialProperty=None, v=None): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ - return self._fastInnerProductDeriv('F', materialProperty=materialProperty) + return self._fastInnerProductDeriv('F', materialProperty=materialProperty, v=v) - def _fastEdgeInnerProductDeriv(self, materialProperty=None): + def _fastEdgeInnerProductDeriv(self, materialProperty=None, v=None): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) :rtype: scipy.csr_matrix :return: M, the inner product matrix (nE, nE) """ - return self._fastInnerProductDeriv('E', materialProperty=materialProperty) + return self._fastInnerProductDeriv('E', materialProperty=materialProperty, v=v) - def _fastInnerProductDeriv(self, AvType, materialProperty=None): + def _fastInnerProductDeriv(self, AvType, materialProperty=None, v=None): """ :param str AvType: 'E' or 'F' :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) @@ -576,11 +576,16 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): """ if materialProperty is None or materialProperty.size == self.nC: Av = getattr(self, 'ave'+AvType+'2CC') - return self.dim * Av.T * Utils.sdiag(self.vol) + V = Utils.sdiag(self.vol) + if v is None: + return self.dim * Av.T * Utils.sdiag(self.vol) + return Utils.sdiag(v) * self.dim * Av.T * V if materialProperty.size == self.nC*self.dim: # anisotropic Av = getattr(self, 'ave'+AvType+'2CCV') V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - return Av.T * V + if v is None: + return Av.T * V + return Utils.sdiag(v) * Av.T * V if __name__ == '__main__': diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py index 80f6bb83..8be5e4eb 100644 --- a/SimPEG/Tests/test_massMatricesDerivs.py +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -5,45 +5,82 @@ from TestUtils import checkDerivative class TestInnerProductsDerivs(unittest.TestCase): - def setUp(self): - pass - def test_FaceIP_derivs_isotropic(self): - for d in range(3): - mesh = Mesh.TensorMesh([10,5,4][d:]) - v = np.random.rand(mesh.nF) - def fun(sig): - M = mesh.getFaceInnerProduct(sig) - Md = mesh.getFaceInnerProductDeriv(sig) - return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(mesh.nC) - passed = checkDerivative(fun, sig, plotIt=False) - self.assertTrue(passed) + def doTestFace(self, h, rep, vec, fast): + mesh = Mesh.TensorMesh(h) + v = np.random.rand(mesh.nF) + def fun(sig): + M = mesh.getFaceInnerProduct(sig) + if vec: + Md = mesh.getFaceInnerProductDeriv(sig, v=v, doFast=fast) + return M*v, Md + Md = mesh.getFaceInnerProductDeriv(sig, doFast=fast) + return M*v, Utils.sdiag(v)*Md + sig = np.random.rand(mesh.nC*rep) + return checkDerivative(fun, sig, num=5, plotIt=False) + + def doTestEdge(self, h, rep, vec, fast): + mesh = Mesh.TensorMesh(h) + v = np.random.rand(mesh.nE) + def fun(sig): + M = mesh.getEdgeInnerProduct(sig) + if vec: + Md = mesh.getEdgeInnerProductDeriv(sig, v=v, doFast=fast) + return M*v, Md + Md = mesh.getEdgeInnerProductDeriv(sig, doFast=fast) + return M*v, Utils.sdiag(v)*Md + sig = np.random.rand(mesh.nC*rep) + return checkDerivative(fun, sig, num=5, plotIt=False) + + def test_FaceIP_1D_isotropic(self): + self.assertTrue(self.doTestFace([10],1,True, False)) + def test_FaceIP_2D_isotropic(self): + self.assertTrue(self.doTestFace([10, 4],1,True, False)) + def test_FaceIP_3D_isotropic(self): + self.assertTrue(self.doTestFace([10, 4, 5],1,True, False)) + def test_FaceIP_2D_anisotropic(self): + self.assertTrue(self.doTestFace([10, 4],2,True, False)) + def test_FaceIP_3D_anisotropic(self): + self.assertTrue(self.doTestFace([10, 4, 5],3,True, False)) + def test_FaceIP_2D_tensor(self): + self.assertTrue(self.doTestFace([10, 4],3,True, False)) + def test_FaceIP_3D_tensor(self): + self.assertTrue(self.doTestFace([10, 4, 5],6,True, False)) + + def test_FaceIP_1D_isotropic_fast(self): + self.assertTrue(self.doTestFace([10],1, False, True)) + def test_FaceIP_2D_isotropic_fast(self): + self.assertTrue(self.doTestFace([10, 4],1, False, True)) + def test_FaceIP_3D_isotropic_fast(self): + self.assertTrue(self.doTestFace([10, 4, 5],1, False, True)) + def test_FaceIP_2D_anisotropic_fast(self): + self.assertTrue(self.doTestFace([10, 4],2, False, True)) + def test_FaceIP_3D_anisotropic_fast(self): + self.assertTrue(self.doTestFace([10, 4, 5],3, False, True)) - def test_EdgeIP_derivs_isotropic(self): - for h in [[10,5],[10,5,4]]: - mesh = Mesh.TensorMesh(h) - v = np.random.rand(mesh.nE) - def fun(sig): - M = mesh.getEdgeInnerProduct(sig) - Md = mesh.getEdgeInnerProductDeriv(sig) - return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(mesh.nC) - passed = checkDerivative(fun, sig, plotIt=False) - self.assertTrue(passed) + def test_EdgeIP_2D_isotropic(self): + self.assertTrue(self.doTestEdge([10, 4],1,True, False)) + def test_EdgeIP_3D_isotropic(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1,True, False)) + def test_EdgeIP_2D_anisotropic(self): + self.assertTrue(self.doTestEdge([10, 4],2,True, False)) + def test_EdgeIP_3D_anisotropic(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3,True, False)) + def test_EdgeIP_2D_tensor(self): + self.assertTrue(self.doTestEdge([10, 4],3,True, False)) + def test_EdgeIP_3D_tensor(self): + self.assertTrue(self.doTestEdge([10, 4, 5],6,True, False)) + + def test_EdgeIP_2D_isotropic_fast(self): + self.assertTrue(self.doTestEdge([10, 4],1, False, True)) + def test_EdgeIP_3D_isotropic_fast(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1, False, True)) + def test_EdgeIP_2D_anisotropic_fast(self): + self.assertTrue(self.doTestEdge([10, 4],2, False, True)) + def test_EdgeIP_3D_anisotropic_fast(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3, False, True)) - def test_FaceIP_derivs_anisotropic(self): - for d in range(3): - mesh = Mesh.TensorMesh([10,5,4][d:]) - v = np.random.rand(mesh.nF) - def fun(sig): - M = mesh.getFaceInnerProduct(sig) - Md = mesh.getFaceInnerProductDeriv(sig) - return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(mesh.nC*mesh.dim) - passed = checkDerivative(fun, sig, plotIt=False) - self.assertTrue(passed) From c6db3d59f75e76dfd8c72ea67439aaa840fdc373 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 3 Mar 2014 11:12:45 -0800 Subject: [PATCH 18/24] refactoring and documentation for innerProducts --- SimPEG/Mesh/InnerProducts.py | 73 +++++++++++++++++++----------------- 1 file changed, 38 insertions(+), 35 deletions(-) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index bb77d5e5..1c85e406 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -74,8 +74,11 @@ class InnerProducts(object): def getFaceInnerProductDeriv(self, materialProperty=None, v=None, P=None, doFast=True): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param numpy.array v: vector to multiply (required in the general implementation) + :param list P: list of projection matrices + :param bool doFast: do a faster implementation if available. :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) + :return: dMdm, the derivative of the inner product matrix (nF, nC*nA) """ fast = None @@ -148,8 +151,11 @@ class InnerProducts(object): def getEdgeInnerProductDeriv(self, materialProperty=None, v=None, P=None, doFast=True): """ :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param numpy.array v: vector to multiply (required in the general implementation) + :param list P: list of projection matrices + :param bool doFast: do a faster implementation if available. :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) + :return: dMdm, the derivative of the inner product matrix (nE, nC*nA) """ fast = None @@ -166,7 +172,14 @@ class InnerProducts(object): return self._getInnerProductDeriv(materialProperty, v, P, self.nE) def _getInnerProductDeriv(self, materialProperty, v, P, n): - + """ + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param numpy.array v: vector to multiply (required in the general implementation) + :param list P: list of projection matrices + :param int n: nF or nE + :rtype: scipy.csr_matrix + :return: dMdm, the derivative of the inner product matrix (n, nC*nA) + """ if v is None: raise Exception('v must be supplied for this implementation.') @@ -186,27 +199,24 @@ class InnerProducts(object): y2 = Y[self.nC:] dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 ))) if materialProperty.size == self.nC*2: - dMdm1 = spzeros(n, self.nC) - dMdm2 = spzeros(n, self.nC) + dMdms = [spzeros(n, self.nC) for _ in range(2)] for i, p in enumerate(P): Y = p * v y1 = Y[:self.nC] y2 = Y[self.nC:] - dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z)) - dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ))) - dMdm = sp.hstack((dMdm1, dMdm2)) + dMdms[0] = dMdms[0] + p.T * sp.vstack(( sdiag( y1 ), Z)) + dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 ))) + dMdm = sp.hstack(dMdms) if materialProperty.size == self.nC*3: - dMdm1 = spzeros(n, self.nC) - dMdm2 = spzeros(n, self.nC) - dMdm3 = spzeros(n, self.nC) + dMdms = [spzeros(n, self.nC) for _ in range(3)] for i, p in enumerate(P): Y = p * v y1 = Y[:self.nC] y2 = Y[self.nC:] - dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z)) - dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ))) - dMdm3 = dMdm3 + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 ))) - dMdm = sp.hstack((dMdm1, dMdm2, dMdm3)) + dMdms[0] = dMdms[0] + p.T * sp.vstack(( sdiag( y1 ), Z)) + dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 ))) + dMdms[2] = dMdms[2] + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 ))) + dMdm = sp.hstack(dMdms) elif d == 3: if materialProperty is None or materialProperty.size == self.nC: dMdm = spzeros(n, self.nC) @@ -217,37 +227,30 @@ class InnerProducts(object): y3 = Y[self.nC*2:] dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 ), sdiag( y3 ))) if materialProperty.size == self.nC*3: - dMdm1 = spzeros(n, self.nC) - dMdm2 = spzeros(n, self.nC) - dMdm3 = spzeros(n, self.nC) + dMdms = [spzeros(n, self.nC) for _ in range(3)] for i, p in enumerate(P): Y = p * v y1 = Y[:self.nC] y2 = Y[self.nC:self.nC*2] y3 = Y[self.nC*2:] - dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z, Z)) - dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ), Z)) - dMdm3 = dMdm3 + p.T * sp.vstack(( Z, Z, sdiag( y3 ))) - dMdm = sp.hstack((dMdm1, dMdm2, dMdm3)) + dMdms[0] = dMdms[0] + p.T * sp.vstack(( sdiag( y1 ), Z, Z)) + dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 ), Z)) + dMdms[2] = dMdms[2] + p.T * sp.vstack(( Z, Z, sdiag( y3 ))) + dMdm = sp.hstack(dMdms) if materialProperty.size == self.nC*6: - dMdm1 = spzeros(n, self.nC) - dMdm2 = spzeros(n, self.nC) - dMdm3 = spzeros(n, self.nC) - dMdm4 = spzeros(n, self.nC) - dMdm5 = spzeros(n, self.nC) - dMdm6 = spzeros(n, self.nC) + dMdms = [spzeros(n, self.nC) for _ in range(6)] for i, p in enumerate(P): Y = p * v y1 = Y[:self.nC] y2 = Y[self.nC:self.nC*2] y3 = Y[self.nC*2:] - dMdm1 = dMdm1 + p.T * sp.vstack(( sdiag( y1 ), Z, Z)) - dMdm2 = dMdm2 + p.T * sp.vstack(( Z, sdiag( y2 ), Z)) - dMdm3 = dMdm3 + p.T * sp.vstack(( Z, Z, sdiag( y3 ))) - dMdm4 = dMdm4 + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 ), Z)) - dMdm5 = dMdm5 + p.T * sp.vstack(( sdiag( y3 ), Z, sdiag( y1 ))) - dMdm6 = dMdm6 + p.T * sp.vstack(( Z, sdiag( y3 ), sdiag( y2 ))) - dMdm = sp.hstack((dMdm1, dMdm2, dMdm3, dMdm4, dMdm5, dMdm6)) + dMdms[0] = dMdms[0] + p.T * sp.vstack(( sdiag( y1 ), Z, Z)) + dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 ), Z)) + dMdms[2] = dMdms[2] + p.T * sp.vstack(( Z, Z, sdiag( y3 ))) + dMdms[3] = dMdms[3] + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 ), Z)) + dMdms[4] = dMdms[4] + p.T * sp.vstack(( sdiag( y3 ), Z, sdiag( y1 ))) + dMdms[5] = dMdms[5] + p.T * sp.vstack(( Z, sdiag( y3 ), sdiag( y2 ))) + dMdm = sp.hstack(dMdms) return dMdm From eeae3ec78365fd7a1f1ad661c100e1108401b0cd Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 3 Mar 2014 12:23:20 -0800 Subject: [PATCH 19/24] 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: From 48e506f4bb12cae57f980a5f89e1ad68212f96e7 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 3 Mar 2014 15:50:06 -0800 Subject: [PATCH 20/24] Clean Plot grid code and add to documentation --- SimPEG/Mesh/LogicallyRectMesh.py | 76 +++++++++++++------------- SimPEG/Mesh/TensorView.py | 94 ++++++++++++-------------------- SimPEG/Mesh/TreeMesh.py | 37 +++++++++---- SimPEG/Utils/meshutils.py | 6 +- docs/api_Mesh.rst | 17 ++++++ docs/api_MeshCode.rst | 18 +++--- 6 files changed, 128 insertions(+), 120 deletions(-) diff --git a/SimPEG/Mesh/LogicallyRectMesh.py b/SimPEG/Mesh/LogicallyRectMesh.py index d606fe78..41dd2641 100644 --- a/SimPEG/Mesh/LogicallyRectMesh.py +++ b/SimPEG/Mesh/LogicallyRectMesh.py @@ -334,7 +334,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts): # Plotting Functions # ############################################# - def plotGrid(self, length=0.05, showIt=False): + def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. @@ -352,55 +352,54 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts): from mpl_toolkits.mplot3d import Axes3D mkvc = Utils.mkvc + axOpts = {'projection':'3d'} if self.dim == 3 else {} + if ax is None: ax = plt.subplot(111, **axOpts) + NN = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: - fig = plt.figure(2) - fig.clf() - ax = plt.subplot(111) - X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten() - X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten() + if lines: + X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten() + Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten() - X = np.r_[X1, X2] - Y = np.r_[Y1, Y2] + X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten() + Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten() - plt.plot(X, Y) + X = np.r_[X1, X2] + Y = np.r_[Y1, Y2] - plt.hold(True) - Nx = self.r(self.normals, 'F', 'Fx', 'V') - Ny = self.r(self.normals, 'F', 'Fy', 'V') - Tx = self.r(self.tangents, 'E', 'Ex', 'V') - Ty = self.r(self.tangents, 'E', 'Ey', 'V') + ax.plot(X, Y, 'b-') + if centers: + ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro') - plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo') + # Nx = self.r(self.normals, 'F', 'Fx', 'V') + # Ny = self.r(self.normals, 'F', 'Fy', 'V') + # Tx = self.r(self.tangents, 'E', 'Ex', 'V') + # Ty = self.r(self.tangents, 'E', 'Ey', 'V') - nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten() - nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten() - plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs') - plt.plot(nX, nY, 'r-') + # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo') - nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten() - nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten() - #plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs') - plt.plot(nX, nY, 'g-') + # nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten() + # nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten() + # ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs') + # ax.plot(nX, nY, 'r-') - tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten() - tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten() - plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^') - plt.plot(tX, tY, 'r-') + # nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten() + # nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten() + # #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs') + # ax.plot(nX, nY, 'g-') - nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten() - nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten() - #plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^') - plt.plot(nX, nY, 'g-') - plt.axis('equal') + # tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten() + # tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten() + # ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^') + # ax.plot(tX, tY, 'r-') + + # nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten() + # nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten() + # #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^') + # ax.plot(nX, nY, 'g-') elif self.dim == 3: - fig = plt.figure(3) - fig.clf() - ax = fig.add_subplot(111, projection='3d') X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten() Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten() Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten() @@ -417,11 +416,10 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts): Y = np.r_[Y1, Y2, Y3] Z = np.r_[Z1, Z2, Z3] - plt.plot(X, Y, 'b', zs=Z) + ax.plot(X, Y, 'b', zs=Z) ax.set_zlabel('x3') ax.grid(True) - ax.hold(False) ax.set_xlabel('x1') ax.set_ylabel('x2') diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index 4c04cda9..3cc98b61 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -369,7 +369,7 @@ class TensorView(object): if showIt: plt.show() return out - def plotGrid(self, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): + def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. :param bool nodes: plot nodes @@ -399,35 +399,26 @@ class TensorView(object): mesh.plotGrid(nodes=True, faces=True, centers=True, lines=True, showIt=True) """ - if self.dim == 1: - fig = plt.figure(1) - fig.clf() - ax = plt.subplot(111) - xn = self.gridN - xc = self.gridCC - ax.hold(True) - ax.plot(xn, np.ones(np.shape(xn)), 'bs') - ax.plot(xc, np.ones(np.shape(xc)), 'ro') - ax.plot(xn, np.ones(np.shape(xn)), 'k--') - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - if showIt: plt.show() - elif self.dim == 2: - fig = plt.figure(2) - fig.clf() - ax = plt.subplot(111) - xn = self.gridN - xc = self.gridCC - xs1 = self.gridFx - xs2 = self.gridFy - ax.hold(True) - if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs') - if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro') + axOpts = {'projection':'3d'} if self.dim == 3 else {} + if ax is None: ax = plt.subplot(111, **axOpts) + + if self.dim == 1: + if nodes: + ax.plot(xn, np.ones(self.nN), 'bs') + if centers: + ax.plot(xc, np.ones(self.nC), 'ro') + if lines: + ax.plot(xn, np.ones(self.nN), 'b-') + ax.set_xlabel('x1') + elif self.dim == 2: + if nodes: + ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bs') + if centers: + ax.plot(self.gridCC[:, 0], self.gridCC[:, 1], 'ro') if faces: - ax.plot(xs1[:, 0], xs1[:, 1], 'g>') - ax.plot(xs2[:, 0], xs2[:, 1], 'g^') + ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'g>') + ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'g^') if edges: ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'c>') ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'c^') @@ -441,38 +432,23 @@ class TensorView(object): Y2 = np.c_[mkvc(NN[1][:, 0]), mkvc(NN[1][:, self.nCy]), mkvc(NN[1][:, 0])*np.nan].flatten() X = np.r_[X1, X2] Y = np.r_[Y1, Y2] - plt.plot(X, Y) + ax.plot(X, Y, 'b-') - ax.grid(True) - ax.hold(False) ax.set_xlabel('x1') ax.set_ylabel('x2') - if showIt: plt.show() elif self.dim == 3: - fig = plt.figure(3) - fig.clf() - ax = fig.add_subplot(111, projection='3d') - xn = self.gridN - xc = self.gridCC - xfs1 = self.gridFx - xfs2 = self.gridFy - xfs3 = self.gridFz - - xes1 = self.gridEx - xes2 = self.gridEy - xes3 = self.gridEz - - ax.hold(True) - if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2]) - if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2]) + if nodes: + ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bs', zs=self.gridN[:, 2]) + if centers: + ax.plot(self.gridCC[:, 0], self.gridCC[:, 1], 'ro', zs=self.gridCC[:, 2]) if faces: - ax.plot(xfs1[:, 0], xfs1[:, 1], 'g>', zs=xfs1[:, 2]) - ax.plot(xfs2[:, 0], xfs2[:, 1], 'g<', zs=xfs2[:, 2]) - ax.plot(xfs3[:, 0], xfs3[:, 1], 'g^', zs=xfs3[:, 2]) + ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'g>', zs=self.gridFx[:, 2]) + ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'g<', zs=self.gridFy[:, 2]) + ax.plot(self.gridFz[:, 0], self.gridFz[:, 1], 'g^', zs=self.gridFz[:, 2]) if edges: - ax.plot(xes1[:, 0], xes1[:, 1], 'k>', zs=xes1[:, 2]) - ax.plot(xes2[:, 0], xes2[:, 1], 'k<', zs=xes2[:, 2]) - ax.plot(xes3[:, 0], xes3[:, 1], 'k^', zs=xes3[:, 2]) + ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'k>', zs=self.gridEx[:, 2]) + ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'k<', zs=self.gridEy[:, 2]) + ax.plot(self.gridEz[:, 0], self.gridEz[:, 1], 'k^', zs=self.gridEz[:, 2]) # Plot the grid lines if lines: @@ -489,14 +465,14 @@ class TensorView(object): X = np.r_[X1, X2, X3] Y = np.r_[Y1, Y2, Y3] Z = np.r_[Z1, Z2, Z3] - plt.plot(X, Y, 'b-', zs=Z) - - ax.grid(True) - ax.hold(False) + ax.plot(X, Y, 'b-', zs=Z) ax.set_xlabel('x1') ax.set_ylabel('x2') ax.set_zlabel('x3') - if showIt: plt.show() + + ax.grid(True) + ax.hold(False) + if showIt: plt.show() def slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None): assert normal in 'xyz', 'normal must be x, y, or z' diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index ae83eae8..8ac841d4 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -322,7 +322,7 @@ class TreeFace(object): if not self.isleaf: return if self.dim == 2: line = np.c_[self.node0.x0, self.node1.x0].T - ax.plot(line[:,0], line[:,1],'r-') + ax.plot(line[:,0], line[:,1],'b-') if text: ax.text(self.center[0], self.center[1],self.num) elif self.dim == 3: if text: ax.text(self.center[0], self.center[1], self.center[2], self.num) @@ -665,10 +665,10 @@ class TreeCell(object): def plotGrid(self, ax, text=False): if not self.isleaf: return if self.dim == 2: - ax.plot(self.center[0],self.center[1],'b.') + ax.plot(self.center[0],self.center[1],'ro') if text: ax.text(self.center[0],self.center[1],self.num) elif self.dim == 3: - ax.plot([self.center[0]],[self.center[1]],'b.', zs=[self.center[2]]) + ax.plot([self.center[0]],[self.center[1]],'ro', zs=[self.center[2]]) if text: ax.text(self.center[0], self.center[1], self.center[2], self.num) @@ -1048,21 +1048,38 @@ class TreeMesh(InnerProducts, BaseMesh): zP = self._getEdgeP(zEdge, xEdge, yEdge) return sp.vstack((xP, yP, zP)) - def plotGrid(self, ax=None, text=True, plotC=True, plotF=True, plotE=False, plotEx=False, plotEy=False, plotEz=False, showIt=False): + def plotGrid(self, ax=None, text=False, centers=False, faces=False, edges=False, lines=True, nodes=False, showIt=False): + self.number() + axOpts = {'projection':'3d'} if self.dim == 3 else {} if ax is None: ax = plt.subplot(111, **axOpts) - if plotC: [c.plotGrid(ax, text=text) for c in self.cells] - if plotF: [f.plotGrid(ax, text=text) for f in self.faces] - if plotE and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edges] - if plotEx and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesX] - if plotEy and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesY] - if plotEz and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesZ] + if lines: + [f.plotGrid(ax, text=text) for f in self.faces] + if centers: + [c.plotGrid(ax, text=text) for c in self.cells] + if faces: + fX = np.array([f.center for f in self.sortedFaceX]) + ax.plot(fX[:,0],fX[:,1],'g>') + fY = np.array([f.center for f in self.sortedFaceY]) + ax.plot(fY[:,0],fY[:,1],'g^') + if edges: + eX = np.array([e.center for e in self.sortedFaceY]) + ax.plot(eX[:,0],eX[:,1],'c>') + eY = np.array([e.center for e in self.sortedFaceX]) + ax.plot(eY[:,0],eY[:,1],'c^') + if nodes: + ns = np.array([n.x0 for n in self.sortedNodes]) + ax.plot(ns[:,0],ns[:,1],'bs') ax.set_xlim((self.x0[0], self.h[0].sum())) ax.set_ylim((self.x0[1], self.h[1].sum())) if self.dim == 3: ax.set_zlim((self.x0[2], self.h[2].sum())) + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') if showIt: plt.show() def plotImage(self, I, ax=None, showIt=True): diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 48dcb8cf..67d747c3 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -11,18 +11,18 @@ def exampleLrmGrid(nC, exType): assert exType in possibleTypes, "Not a possible example type." if exType == 'rect': - return ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) + return list(ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)) elif exType == 'rotate': if len(nC) == 2: X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2) amt[amt < 0] = 0 - return X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt + return [X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt] elif len(nC) == 3: X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False) amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2) amt[amt < 0] = 0 - return X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt + return [X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt] def meshTensors(*args): diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst index 46bd877b..e5a72491 100644 --- a/docs/api_Mesh.rst +++ b/docs/api_Mesh.rst @@ -21,7 +21,24 @@ Each mesh code follows the guiding principles that are present in this tutorial, but the details, advantages and disadvantages differ between the implementations. +.. plot:: + from SimPEG import Mesh, Utils, np + import matplotlib.pyplot as plt + sz = [10,10] + tM = Mesh.TensorMesh(sz) + qM = Mesh.TreeMesh(sz) + qM.refine(lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.3 else 0) + rM = Mesh.LogicallyRectMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) + + fig, axes = plt.subplots(1,3,figsize=(14,5)) + opts = {} + tM.plotGrid(ax=axes[0], **opts) + axes[0].set_title('TensorMesh') + qM.plotGrid(ax=axes[1], **opts) + axes[1].set_title('TreeMesh') + rM.plotGrid(ax=axes[2], **opts) + axes[2].set_title('LogicallyRectMesh') Variable Locations and Terminology diff --git a/docs/api_MeshCode.rst b/docs/api_MeshCode.rst index 090bb6c2..65820fc8 100644 --- a/docs/api_MeshCode.rst +++ b/docs/api_MeshCode.rst @@ -9,15 +9,6 @@ Tensor Mesh :undoc-members: -Cylindrical 1D Mesh -=================== - -.. automodule:: SimPEG.Mesh.Cyl1DMesh - :show-inheritance: - :members: - :undoc-members: - - Logically Rectangular Mesh ========================== @@ -27,6 +18,15 @@ Logically Rectangular Mesh :undoc-members: +Cylindrical 1D Mesh +=================== + +.. automodule:: SimPEG.Mesh.Cyl1DMesh + :show-inheritance: + :members: + :undoc-members: + + Base Mesh ========= From e59a0390522ff8698d51883b568200c2e9bf8a77 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 4 Mar 2014 09:56:12 -0800 Subject: [PATCH 21/24] Float issue on inner products. --- SimPEG/Mesh/TensorMesh.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index e142d78d..df049970 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -534,6 +534,8 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): """ if materialProperty is None: materialProperty = np.ones(self.nC) + elif type(materialProperty) in [float, int, long]: + materialProperty = materialProperty * np.ones(M.nC) if materialProperty.size == self.nC: if invertProperty: From 90967ce7f976b993565ff56a3ab8e2126ffc6e50 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 4 Mar 2014 12:19:57 -0800 Subject: [PATCH 22/24] isScalar in Utils (tested) and support for scalar inner products and derivatives (tested). --- SimPEG/Mesh/InnerProducts.py | 28 +++++++++++++++--------- SimPEG/Mesh/TensorMesh.py | 29 ++++++++++++++++--------- SimPEG/Tests/test_massMatrices.py | 2 +- SimPEG/Tests/test_massMatricesDerivs.py | 24 ++++++++++++++++++-- SimPEG/Tests/test_utils.py | 6 +++++ SimPEG/Utils/matutils.py | 13 +++++++++-- 6 files changed, 77 insertions(+), 25 deletions(-) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 99490674..22c6a521 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -1,5 +1,5 @@ from scipy import sparse as sp -from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor, invPropertyTensor, spzeros +from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor, invPropertyTensor, spzeros, isScalar import numpy as np @@ -183,25 +183,33 @@ class InnerProducts(object): :rtype: scipy.csr_matrix :return: dMdm, the derivative of the inner product matrix (n, nC*nA) """ + if materialProperty is None: + return None + if v is None: raise Exception('v must be supplied for this implementation.') d = self.dim Z = spzeros(self.nC, self.nC) - if d == 1: - dMdm = spzeros(n, self.nC) + if isScalar(materialProperty): + dMdm = spzeros(n, 1) for i, p in enumerate(P): - dMdm = dMdm + p.T * sdiag( p * v ) + dMdm = dMdm + sp.csr_matrix((p.T * (p * v), (range(n), np.zeros(n))), shape=(n,1)) + if d == 1: + if materialProperty.size == self.nC: + dMdm = spzeros(n, self.nC) + for i, p in enumerate(P): + dMdm = dMdm + p.T * sdiag( p * v ) elif d == 2: - if materialProperty is None or materialProperty.size == self.nC: + if materialProperty.size == self.nC: dMdm = spzeros(n, self.nC) for i, p in enumerate(P): Y = p * v y1 = Y[:self.nC] y2 = Y[self.nC:] dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 ))) - if materialProperty.size == self.nC*2: + elif materialProperty.size == self.nC*2: dMdms = [spzeros(n, self.nC) for _ in range(2)] for i, p in enumerate(P): Y = p * v @@ -210,7 +218,7 @@ class InnerProducts(object): dMdms[0] = dMdms[0] + p.T * sp.vstack(( sdiag( y1 ), Z)) dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 ))) dMdm = sp.hstack(dMdms) - if materialProperty.size == self.nC*3: + elif materialProperty.size == self.nC*3: dMdms = [spzeros(n, self.nC) for _ in range(3)] for i, p in enumerate(P): Y = p * v @@ -221,7 +229,7 @@ class InnerProducts(object): dMdms[2] = dMdms[2] + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 ))) dMdm = sp.hstack(dMdms) elif d == 3: - if materialProperty is None or materialProperty.size == self.nC: + if materialProperty.size == self.nC: dMdm = spzeros(n, self.nC) for i, p in enumerate(P): Y = p * v @@ -229,7 +237,7 @@ class InnerProducts(object): y2 = Y[self.nC:self.nC*2] y3 = Y[self.nC*2:] dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 ), sdiag( y3 ))) - if materialProperty.size == self.nC*3: + elif materialProperty.size == self.nC*3: dMdms = [spzeros(n, self.nC) for _ in range(3)] for i, p in enumerate(P): Y = p * v @@ -240,7 +248,7 @@ class InnerProducts(object): dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 ), Z)) dMdms[2] = dMdms[2] + p.T * sp.vstack(( Z, Z, sdiag( y3 ))) dMdm = sp.hstack(dMdms) - if materialProperty.size == self.nC*6: + elif materialProperty.size == self.nC*6: dMdms = [spzeros(n, self.nC) for _ in range(6)] for i, p in enumerate(P): Y = p * v diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index df049970..1bc7cadb 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -534,18 +534,18 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): """ if materialProperty is None: materialProperty = np.ones(self.nC) - elif type(materialProperty) in [float, int, long]: - materialProperty = materialProperty * np.ones(M.nC) + + if invertProperty: + materialProperty = 1./materialProperty + + if Utils.isScalar(materialProperty): + materialProperty = materialProperty*np.ones(self.nC) if materialProperty.size == self.nC: - if invertProperty: - materialProperty = 1./materialProperty Av = getattr(self, 'ave'+AvType+'2CC') - V = Utils.sdiag(self.vol) - return self.dim * Utils.sdiag(Av.T * V * materialProperty) + Vprop = self.vol * Utils.mkvc(materialProperty) + return self.dim * Utils.sdiag(Av.T * Vprop) if materialProperty.size == self.nC*self.dim: - if invertProperty: - materialProperty = 1./materialProperty Av = getattr(self, 'ave'+AvType+'2CCV') V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty)) @@ -576,11 +576,20 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): :rtype: scipy.csr_matrix :return: M, the inner product matrix (nF, nF) """ - if materialProperty is None or materialProperty.size == self.nC: + if materialProperty is None: + return None + if Utils.isScalar(materialProperty): + Av = getattr(self, 'ave'+AvType+'2CC') + V = Utils.sdiag(self.vol) + ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1)) + if v is None: + return self.dim * Av.T * V * ones + return Utils.sdiag(v) * self.dim * Av.T * V * ones + if materialProperty.size == self.nC: Av = getattr(self, 'ave'+AvType+'2CC') V = Utils.sdiag(self.vol) if v is None: - return self.dim * Av.T * Utils.sdiag(self.vol) + return self.dim * Av.T * V return Utils.sdiag(v) * self.dim * Av.T * V if materialProperty.size == self.nC*self.dim: # anisotropic Av = getattr(self, 'ave'+AvType+'2CCV') diff --git a/SimPEG/Tests/test_massMatrices.py b/SimPEG/Tests/test_massMatrices.py index 3d3946ad..67de753e 100644 --- a/SimPEG/Tests/test_massMatrices.py +++ b/SimPEG/Tests/test_massMatrices.py @@ -122,7 +122,7 @@ class TestInnerProducts2D(OrderTest): sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)] analytic = 189959./120 # Found using sympy. z=5 elif self.sigmaTest == 3: - sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] + sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] analytic = 781427./360 # Found using sympy. z=5 if self.location == 'edges': diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py index 8be5e4eb..d66c0b65 100644 --- a/SimPEG/Tests/test_massMatricesDerivs.py +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -16,7 +16,7 @@ class TestInnerProductsDerivs(unittest.TestCase): return M*v, Md Md = mesh.getFaceInnerProductDeriv(sig, doFast=fast) return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(mesh.nC*rep) + sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep) return checkDerivative(fun, sig, num=5, plotIt=False) def doTestEdge(self, h, rep, vec, fast): @@ -29,9 +29,15 @@ class TestInnerProductsDerivs(unittest.TestCase): return M*v, Md Md = mesh.getEdgeInnerProductDeriv(sig, doFast=fast) return M*v, Utils.sdiag(v)*Md - sig = np.random.rand(mesh.nC*rep) + sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep) return checkDerivative(fun, sig, num=5, plotIt=False) + def test_FaceIP_1D_float(self): + self.assertTrue(self.doTestFace([10],0,True, False)) + def test_FaceIP_2D_float(self): + self.assertTrue(self.doTestFace([10, 4],0,True, False)) + def test_FaceIP_3D_float(self): + self.assertTrue(self.doTestFace([10, 4, 5],0,True, False)) def test_FaceIP_1D_isotropic(self): self.assertTrue(self.doTestFace([10],1,True, False)) def test_FaceIP_2D_isotropic(self): @@ -47,6 +53,12 @@ class TestInnerProductsDerivs(unittest.TestCase): def test_FaceIP_3D_tensor(self): self.assertTrue(self.doTestFace([10, 4, 5],6,True, False)) + def test_FaceIP_1D_float_fast(self): + self.assertTrue(self.doTestFace([10],0, False, True)) + def test_FaceIP_2D_float_fast(self): + self.assertTrue(self.doTestFace([10, 4],0, False, True)) + def test_FaceIP_3D_float_fast(self): + self.assertTrue(self.doTestFace([10, 4, 5],0, False, True)) def test_FaceIP_1D_isotropic_fast(self): self.assertTrue(self.doTestFace([10],1, False, True)) def test_FaceIP_2D_isotropic_fast(self): @@ -59,6 +71,10 @@ class TestInnerProductsDerivs(unittest.TestCase): self.assertTrue(self.doTestFace([10, 4, 5],3, False, True)) + def test_EdgeIP_2D_float(self): + self.assertTrue(self.doTestEdge([10, 4],0,True, False)) + def test_EdgeIP_3D_float(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0,True, False)) def test_EdgeIP_2D_isotropic(self): self.assertTrue(self.doTestEdge([10, 4],1,True, False)) def test_EdgeIP_3D_isotropic(self): @@ -72,6 +88,10 @@ class TestInnerProductsDerivs(unittest.TestCase): def test_EdgeIP_3D_tensor(self): self.assertTrue(self.doTestEdge([10, 4, 5],6,True, False)) + def test_EdgeIP_2D_float_fast(self): + self.assertTrue(self.doTestEdge([10, 4],0, False, True)) + def test_EdgeIP_3D_float_fast(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0, False, True)) def test_EdgeIP_2D_isotropic_fast(self): self.assertTrue(self.doTestEdge([10, 4],1, False, True)) def test_EdgeIP_3D_isotropic_fast(self): diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index e69c8473..ad737489 100644 --- a/SimPEG/Tests/test_utils.py +++ b/SimPEG/Tests/test_utils.py @@ -164,6 +164,12 @@ class TestSequenceFunctions(unittest.TestCase): Z = B2*A - sp.identity(M.nC*3) self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL) + def test_isFloat(self): + self.assertTrue(isScalar(1.)) + self.assertTrue(isScalar(1)) + self.assertTrue(isScalar(long(1))) + self.assertTrue(isScalar(np.r_[1.])) + self.assertTrue(isScalar(np.r_[1])) if __name__ == '__main__': unittest.main() diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index c15f6234..f1ea7ec7 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -1,6 +1,15 @@ import numpy as np import scipy.sparse as sp + +def isScalar(f): + scalarTypes = [float, int, long, np.float_, np.int_] + if type(f) in scalarTypes: + return True + elif type(f) == np.ndarray and f.size == 1 and type(f[0]) in scalarTypes: + return True + return False + def mkvc(x, numDims=1): """Creates a vector with the number of dimension specified @@ -248,7 +257,7 @@ def makePropertyTensor(M, sigma): if sigma is None: # default is ones sigma = np.ones(M.nC) - if type(sigma) in [float, int, long]: + if isScalar(sigma): sigma = sigma * np.ones(M.nC) if M.dim == 1: @@ -294,7 +303,7 @@ def invPropertyTensor(M, tensor, returnMatrix=False): T = None - if type(tensor) in [float, int, long]: + if isScalar(tensor): T = 1./tensor elif tensor.size == M.nC: # Isotropic! From 95bbfb8d13178cf673d8ac4b8a586798308c498e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 5 Mar 2014 16:20:55 -0800 Subject: [PATCH 23/24] Updates to check derivative code (should now work for complex, without throwing warnings.) --- SimPEG/Tests/TestUtils.py | 37 +++++++++++++++++++++++-------------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index da18f537..c4d11443 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -220,35 +220,44 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole """ print "%s checkDerivative %s" % ('='*20, '='*20) - print "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|\tOrder\n%s" % ('-'*57) + print "iter h |f0-ft| |f0-ft-h*J0*dx| Order\n%s" % ('-'*57) - Jc = fctn(x0) + f0, J0 = fctn(x0) x0 = mkvc(x0) if dx is None: dx = np.random.randn(len(x0)) - t = np.logspace(-1, -num, num) - E0 = np.ones(t.shape) - E1 = np.ones(t.shape) + h = np.logspace(-1, -num, num) + E0 = np.ones(h.shape) + E1 = np.ones(h.shape) + + def l2norm(x): + # because np.norm breaks if they are scalars? + return np.sqrt(np.real(np.vdot(x, x))) - l2norm = lambda x: np.sqrt(np.inner(x, x)) # because np.norm breaks if they are scalars? for i in range(num): - Jt = fctn(x0+t[i]*dx) - E0[i] = l2norm(Jt[0]-Jc[0]) # 0th order Taylor - if inspect.isfunction(Jc[1]): - E1[i] = l2norm(Jt[0]-Jc[0]-t[i]*Jc[1](dx)) # 1st order Taylor + # Evaluate at test point + ft, Jt = fctn( x0 + h[i]*dx ) + # 0th order Taylor + E0[i] = l2norm( ft - f0 ) + # 1st order Taylor + if inspect.isfunction(J0): + E1[i] = l2norm( ft - f0 - h[i]*J0(dx) ) else: # We assume it is a numpy.ndarray - E1[i] = l2norm(Jt[0]-Jc[0]-t[i]*Jc[1].dot(dx)) # 1st order Taylor + E1[i] = l2norm( ft - f0 - h[i]*J0.dot(dx) ) + order0 = np.log10(E0[:-1]/E0[1:]) order1 = np.log10(E1[:-1]/E1[1:]) - print "%d\t%1.2e\t%1.3e\t\t%1.3e\t\t%1.3f" % (i, t[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1]) + print " %d %1.2e %1.3e %1.3e %1.3f" % (i, h[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1]) + # Ensure we are about precision order0 = order0[E0[1:] > eps] order1 = order1[E1[1:] > eps] belowTol = order1.size == 0 and order0.size > 0 + # Make sure we get the correct order correctOrder = order1.size > 0 and np.mean(order1) > tolerance * expectedOrder passTest = belowTol or correctOrder @@ -264,8 +273,8 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole if plotIt: plt.figure() plt.clf() - plt.loglog(t, E0, 'b') - plt.loglog(t, E1, 'g--') + plt.loglog(h, E0, 'b') + plt.loglog(h, E1, 'g--') plt.title('checkDerivative') plt.xlabel('h') plt.ylabel('error of Taylor approximation') From a6badc51c972cef7aa4adaffd571154d7bd82ce6 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 5 Mar 2014 16:21:27 -0800 Subject: [PATCH 24/24] bug fix for tensorview plotSlice --- SimPEG/Mesh/TensorView.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index 3cc98b61..a67d4609 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -305,7 +305,7 @@ class TensorView(object): # Now just deal with 'F' and 'E' aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC') v = getattr(self,aveOp)*v # average to cell centers (might be a vector) - v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M') + v = self.r(v.reshape((self.nC,-1),order='F'),'CC','CC','M') if view == 'vec': outSlice = [] if 'X' not in normal: outSlice.append(getIndSlice(v[0]))