mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-10 07:48:48 +08:00
inner product documentation
This commit is contained in:
Rowan Cockett
+193
-4
@@ -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))
|
||||
|
||||
@@ -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):
|
||||
|
||||
@@ -1,8 +1,8 @@
|
||||
.. _api_LOMView:
|
||||
|
||||
LOM View
|
||||
***********
|
||||
********
|
||||
|
||||
.. automodule:: SimPEG.LOMView
|
||||
.. automodule:: SimPEG.LomView
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
+1
-1
@@ -26,7 +26,7 @@ Meshing & Operators
|
||||
api_TensorMesh
|
||||
api_TensorView
|
||||
api_LogicallyOrthogonalMesh
|
||||
api_LomView
|
||||
api_LOMView
|
||||
api_DiffOperators
|
||||
api_InnerProducts
|
||||
|
||||
|
||||
Reference in New Issue
Block a user