diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py index 1add3e72..18803344 100644 --- a/SimPEG/InnerProducts.py +++ b/SimPEG/InnerProducts.py @@ -43,6 +43,88 @@ class InnerProducts(object): def getFaceInnerProduct(mesh, mu=None, returnP=False): + """ + :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :rtype: scipy.csr_matrix + :return: M, the inner product matrix + + Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows: + + .. math:: + \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{1} & 0 \\\\ 0 & 0 & \mu_{1} \end{matrix}\\right] + + \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{2} & 0 \\\\ 0 & 0 & \mu_{3} \end{matrix}\\right] + + \left[\\begin{matrix} \mu_{1} & \mu_{4} & \mu_{5} \\\\ \mu_{4} & \mu_{2} & \mu_{6} \\\\ \mu_{5} & \mu_{6} & \mu_{3} \end{matrix}\\right] + + Example problem for DC resistivity: + + .. math:: + + \sigma^{-1}\mathbf{J} = \\nabla \phi + + We can define in weak form by integrating with a general face function F: + + .. math:: + + \int_{\\text{cell}}{\sigma^{-1}\mathbf{J} \cdot \mathbf{F}} = \int_{\\text{cell}}{\\nabla \phi \cdot \mathbf{F}} + + \int_{\\text{cell}}{\sigma^{-1}\mathbf{J} \cdot \mathbf{F}} = \int_{\\text{cell}}{(\\nabla \cdot \mathbf{F}) \phi } + \int_{\partial \\text{cell}}{ \phi \mathbf{F} \cdot \mathbf{n}} + + We can then discretize for every cell: + + .. math:: + + v_{\\text{cell}} \sigma^{-1} (\mathbf{J}_x \mathbf{F}_x +\mathbf{J}_y \mathbf{F}_y + \mathbf{J}_z \mathbf{F}_z ) = -\phi^{\\top} v_{\\text{cell}} (\mathbf{D}_{\\text{cell}} \mathbf{F}) + \\text{BC} + + We can represent this in vector form (again this is for every cell), and will generalize for the case of anisotropic (tensor) sigma. + + .. math:: + + \mathbf{F}_c^{\\top} (\sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}}) \mathbf{J}_c = -\phi^{\\top} v_{\\text{cell}}( v_\\text{cell}^{-1} \mathbf{D}_{\\text{cell}} \mathbf{A} \mathbf{F}) + \\text{BC} + + We multiply by volume on each side of the tensor conductivity to keep symmetry in the system. Here J_c is the Cartesian J (on the faces) and must be calculated differently depending on the mesh: + + .. math:: + \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\\text{TENSOR} = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\\text{LOM} + + Here the i index refers to where we choose to approximate this integral. We will approximate this relation at every node of the cell, there are 8 in 3D, using a projection matrix Q_i to pick the appropriate fluxes. We will then average to the cell center: + + .. math:: + + \mathbf{F}^{\\top} + {1\over 8} + \left(\sum_{i=1}^8 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + \mathbf{J} + = + -\mathbf{F}^{\\top} \mathbf{A} \mathbf{D}_{\\text{cell}}^{\\top} \phi + \\text{BC} + + \mathbf{M}(\Sigma^{-1}) \mathbf{J} + = + -\mathbf{A} \mathbf{D}_{\\text{cell}}^{\\top} \phi + \\text{BC} + + \mathbf{M}(\Sigma^{-1}) = {1\over 8} + \left(\sum_{i=1}^8 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + The M is returned if mu is set equal to \Sigma^{-1}. + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P000, P001, P010, P011, P100, P101, P110, P111] + + Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + + """ if mu is None: # default is ones mu = np.ones((mesh.nC, 1)) @@ -82,10 +164,10 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False): # 100 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k # 010 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k # 110 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k - # 001 | i ,j ,k | i , j, k | i, j , k | i, j, k+1 - # 101 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k+1 - # 011 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k+1 - # 111 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k+1 + # 001 | i ,j ,k+1 | i , j, k | i, j , k | i, j, k+1 + # 101 | i+1,j ,k+1 | i+1, j, k | i, j , k | i, j, k+1 + # 011 | i ,j+1,k+1 | i , j, k | i, j+1, k | i, j, k+1 + # 111 | i+1,j+1,k+1 | i+1, j, k | i, j+1, k | i, j, k+1 # Square root of cell volume multiplied by 1/8 v = np.sqrt(0.125*mesh.vol) @@ -120,6 +202,42 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False): def getFaceInnerProduct2D(mesh, mu=None, returnP=False): + """ + :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 2, or 3)) + :param bool returnP: returns the projection matrices + :rtype: scipy.csr_matrix + :return: M, the inner product matrix + + Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows: + + .. math:: + \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{1} \end{matrix}\\right] + + \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{2} \end{matrix}\\right] + + \left[\\begin{matrix} \mu_{1} & \mu_{3} \\\\ \mu_{3} & \mu_{2} \end{matrix}\\right] + + + .. math:: + + \mathbf{M}(\Sigma^{-1}) = {1\over 4} + \left(\sum_{i=1}^4 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P00, P10, P01, P11] + + Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + + """ if mu is None: # default is ones mu = np.ones((mesh.nC, 1)) @@ -185,6 +303,41 @@ def getFaceInnerProduct2D(mesh, mu=None, returnP=False): def getEdgeInnerProduct(mesh, sigma=None, returnP=False): + """ + :param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :rtype: scipy.csr_matrix + :return: M, the inner product matrix + + + Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows: + + .. math:: + \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{1} & 0 \\\\ 0 & 0 & \sigma_{1} \end{matrix}\\right] + + \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{2} & 0 \\\\ 0 & 0 & \sigma_{3} \end{matrix}\\right] + + \left[\\begin{matrix} \sigma_{1} & \sigma_{4} & \sigma_{5} \\\\ \sigma_{4} & \sigma_{2} & \sigma_{6} \\\\ \sigma_{5} & \sigma_{6} & \sigma_{3} \end{matrix}\\right] + + What is returned: + + .. math:: + \mathbf{M}(\Sigma^{-1}) = {1\over 8} + \left(\sum_{i=1}^8 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P000, P001, P010, P011, P100, P101, P110, P111] + + Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + """ if sigma is None: # default is ones sigma = np.ones((mesh.nC, 1)) @@ -262,6 +415,42 @@ def getEdgeInnerProduct(mesh, sigma=None, returnP=False): def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False): + """ + :param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 2, or 3)) + :param bool returnP: returns the projection matrices + :rtype: scipy.csr_matrix + :return: M, the inner product matrix + + Depending on the number of columns (either 1, 2, or 3) of sigma, the material property is interpreted as follows: + + .. math:: + \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{1} \end{matrix}\\right] + + \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{2} \end{matrix}\\right] + + \left[\\begin{matrix} \sigma_{1} & \sigma_{3} \\\\ \sigma_{3} & \sigma_{2} \end{matrix}\\right] + + + .. math:: + + \mathbf{M}(\Sigma^{-1}) = {1\over 4} + \left(\sum_{i=1}^4 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P00, P10, P01, P11] + + Here each P is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + + """ if sigma is None: # default is ones sigma = np.ones((mesh.nC, 1)) diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/LogicallyOrthogonalMesh.py index 49e2ee07..b9516a6d 100644 --- a/SimPEG/LogicallyOrthogonalMesh.py +++ b/SimPEG/LogicallyOrthogonalMesh.py @@ -258,17 +258,16 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): area = property(**area()) def normals(): - doc = """ -Face normals: calling this will average -the computed normals so that there is one -per face. This is especially relevant in -3D, as there are up to 4 different normals -for each face that will be different. + doc = """Face normals: calling this will average + the computed normals so that there is one + per face. This is especially relevant in + 3D, as there are up to 4 different normals + for each face that will be different. -To reshape the normals into a matrix and get the y component: + To reshape the normals into a matrix and get the y component: -NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M') -""" + NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M') + """ def fget(self): if(self._normals is None): diff --git a/docs/api_LOMView.rst b/docs/api_LOMView.rst index 67e1bb91..61630c26 100644 --- a/docs/api_LOMView.rst +++ b/docs/api_LOMView.rst @@ -1,8 +1,8 @@ .. _api_LOMView: LOM View -*********** +******** -.. automodule:: SimPEG.LOMView +.. automodule:: SimPEG.LomView :members: :undoc-members: diff --git a/docs/index.rst b/docs/index.rst index ce50bc3d..22c9d5e9 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -26,7 +26,7 @@ Meshing & Operators api_TensorMesh api_TensorView api_LogicallyOrthogonalMesh - api_LomView + api_LOMView api_DiffOperators api_InnerProducts