mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-28 11:26:12 +08:00
Merge branch 'master' of https://github.com/simpeg/simpeg into cylClean
Conflicts: SimPEG/Mesh/LogicallyRectMesh.py SimPEG/Mesh/TensorMesh.py SimPEG/Mesh/__init__.py SimPEG/Tests/TestUtils.py SimPEG/Tests/test_operators.py
This commit is contained in:
rowanc1
@@ -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):
|
||||
|
||||
@@ -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
|
||||
@@ -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
|
||||
|
||||
+225
-254
@@ -1,187 +1,70 @@
|
||||
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, spzeros, isScalar
|
||||
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,
|
||||
invertProperty=False, doFast=True):
|
||||
"""
|
||||
: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
|
||||
: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 (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
|
||||
fast = None
|
||||
|
||||
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 returnP is False and hasattr(self, '_fastFaceInnerProduct') and doFast:
|
||||
fast = self._fastFaceInnerProduct(materialProperty=materialProperty, invertProperty=invertProperty)
|
||||
|
||||
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
|
||||
if fast is not None:
|
||||
return fast
|
||||
|
||||
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')
|
||||
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)))
|
||||
|
||||
if d == 1:
|
||||
fP = _getFacePx(self)
|
||||
P000 = V*fP('fXm')
|
||||
P100 = V*fP('fXp')
|
||||
elif d == 2:
|
||||
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:
|
||||
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(M, mu)
|
||||
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 +72,65 @@ class InnerProducts(object):
|
||||
else:
|
||||
return A
|
||||
|
||||
def getEdgeInnerProduct(M, sigma=None, returnP=False):
|
||||
def getFaceInnerProductDeriv(self, materialProperty=None, v=None, P=None, doFast=True):
|
||||
"""
|
||||
: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
|
||||
: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 (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: dMdm, the derivative of the inner product matrix (nF, nC*nA)
|
||||
"""
|
||||
if M.dim == 1:
|
||||
fast = None
|
||||
|
||||
if hasattr(self, '_fastFaceInnerProductDeriv') and doFast:
|
||||
fast = self._fastFaceInnerProductDeriv(materialProperty=materialProperty, v=v)
|
||||
|
||||
if fast is not None:
|
||||
return fast
|
||||
|
||||
if P is None:
|
||||
M, P = self.getFaceInnerProduct(materialProperty=materialProperty, returnP=True)
|
||||
|
||||
return self._getInnerProductDeriv(materialProperty, v, P, self.nF)
|
||||
|
||||
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)
|
||||
"""
|
||||
fast = None
|
||||
|
||||
if returnP is False and hasattr(self, '_fastEdgeInnerProduct') and doFast:
|
||||
fast = self._fastEdgeInnerProduct(materialProperty=materialProperty, invertProperty=invertProperty)
|
||||
|
||||
if fast is not None:
|
||||
return fast
|
||||
|
||||
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)))
|
||||
|
||||
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,17 +140,131 @@ 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
|
||||
else:
|
||||
return A
|
||||
|
||||
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: dMdm, the derivative of the inner product matrix (nE, nC*nA)
|
||||
"""
|
||||
|
||||
fast = None
|
||||
|
||||
if hasattr(self, '_fastEdgeInnerProductDeriv') and doFast:
|
||||
fast = self._fastEdgeInnerProductDeriv(materialProperty=materialProperty, v=v)
|
||||
|
||||
if fast is not None:
|
||||
return fast
|
||||
|
||||
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):
|
||||
"""
|
||||
: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 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 isScalar(materialProperty):
|
||||
dMdm = spzeros(n, 1)
|
||||
for i, p in enumerate(P):
|
||||
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.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 )))
|
||||
elif materialProperty.size == self.nC*2:
|
||||
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:]
|
||||
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)
|
||||
elif materialProperty.size == self.nC*3:
|
||||
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:]
|
||||
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.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 )))
|
||||
elif materialProperty.size == self.nC*3:
|
||||
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:]
|
||||
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)
|
||||
elif materialProperty.size == self.nC*6:
|
||||
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:]
|
||||
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
|
||||
|
||||
# ------------------------ Geometries ------------------------------
|
||||
#
|
||||
#
|
||||
@@ -380,11 +351,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))
|
||||
@@ -392,7 +363,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')
|
||||
|
||||
@@ -417,7 +388,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
|
||||
@@ -440,7 +411,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')
|
||||
@@ -474,7 +445,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]))
|
||||
@@ -489,7 +460,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')
|
||||
|
||||
@@ -509,7 +480,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
|
||||
@@ -523,7 +494,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')
|
||||
@@ -552,7 +523,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]))
|
||||
|
||||
@@ -2,7 +2,6 @@ from SimPEG import Utils, np
|
||||
from BaseMesh import BaseRectangularMesh
|
||||
from DiffOperators import DiffOperators
|
||||
from InnerProducts import InnerProducts
|
||||
from View 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,104 @@ class LogicallyOrthogonalMesh(BaseRectangularMesh, DiffOperators, InnerProducts,
|
||||
_tangents = None
|
||||
tangents = property(**tangents())
|
||||
|
||||
|
||||
|
||||
#############################################
|
||||
# Plotting Functions #
|
||||
#############################################
|
||||
|
||||
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.
|
||||
|
||||
|
||||
.. 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
|
||||
|
||||
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:
|
||||
|
||||
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()
|
||||
|
||||
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]
|
||||
|
||||
ax.plot(X, Y, 'b-')
|
||||
if centers:
|
||||
ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
|
||||
|
||||
# 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(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()
|
||||
# ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
|
||||
# ax.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()
|
||||
# #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
|
||||
# ax.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()
|
||||
# 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:
|
||||
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]
|
||||
|
||||
ax.plot(X, Y, 'b', zs=Z)
|
||||
ax.set_zlabel('x3')
|
||||
|
||||
ax.grid(True)
|
||||
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 +435,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')
|
||||
@@ -20,6 +20,7 @@ class BaseTensorMesh(BaseRectangularMesh):
|
||||
|
||||
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, np.int_]:
|
||||
@@ -445,6 +446,112 @@ class TensorMesh(BaseTensorMesh, 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)
|
||||
"""
|
||||
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 invertProperty:
|
||||
materialProperty = 1./materialProperty
|
||||
|
||||
if Utils.isScalar(materialProperty):
|
||||
materialProperty = materialProperty*np.ones(self.nC)
|
||||
|
||||
if materialProperty.size == self.nC:
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
Vprop = self.vol * Utils.mkvc(materialProperty)
|
||||
return self.dim * Utils.sdiag(Av.T * Vprop)
|
||||
if materialProperty.size == self.nC*self.dim:
|
||||
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, 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, v=v)
|
||||
|
||||
|
||||
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, v=v)
|
||||
|
||||
|
||||
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))
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
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 * V
|
||||
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))
|
||||
if v is None:
|
||||
return Av.T * V
|
||||
return Utils.sdiag(v) * Av.T * V
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
print('Welcome to tensor mesh!')
|
||||
|
||||
+27
-10
@@ -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):
|
||||
|
||||
+36
-60
@@ -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]))
|
||||
@@ -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'
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
from TensorMesh import TensorMesh
|
||||
from CylMesh import CylMesh
|
||||
from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
|
||||
from LogicallyRectMesh import LogicallyRectMesh
|
||||
from TreeMesh import TreeMesh
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
+10
-4
@@ -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'
|
||||
|
||||
@@ -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)
|
||||
|
||||
|
||||
|
||||
+3
-1
@@ -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
|
||||
|
||||
@@ -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) ) )
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
+1
-1
@@ -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:
|
||||
|
||||
+29
-20
@@ -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, CylMesh
|
||||
from SimPEG.Mesh import TensorMesh, LogicallyRectMesh, CylMesh
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import unittest
|
||||
@@ -115,7 +115,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:
|
||||
@@ -125,11 +125,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):
|
||||
@@ -231,35 +231,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
|
||||
@@ -275,8 +284,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')
|
||||
|
||||
@@ -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()
|
||||
@@ -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()
|
||||
@@ -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__':
|
||||
|
||||
@@ -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]
|
||||
|
||||
@@ -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),
|
||||
@@ -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]
|
||||
|
||||
@@ -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':
|
||||
|
||||
@@ -0,0 +1,108 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from TestUtils import checkDerivative
|
||||
|
||||
|
||||
class TestInnerProductsDerivs(unittest.TestCase):
|
||||
|
||||
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(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):
|
||||
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(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):
|
||||
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_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):
|
||||
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_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):
|
||||
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_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):
|
||||
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))
|
||||
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -4,7 +4,7 @@ from TestUtils import OrderTest
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
#TODO: 'randomTensorMesh'
|
||||
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)]
|
||||
@@ -38,7 +38,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)
|
||||
@@ -208,7 +208,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)
|
||||
|
||||
@@ -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):
|
||||
@@ -163,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()
|
||||
|
||||
@@ -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
|
||||
|
||||
+31
-13
@@ -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
|
||||
|
||||
@@ -146,7 +155,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')
|
||||
@@ -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:
|
||||
@@ -261,9 +270,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 +284,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
|
||||
|
||||
|
||||
@@ -289,34 +303,38 @@ 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!
|
||||
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
|
||||
|
||||
@@ -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()
|
||||
@@ -11,18 +11,18 @@ def exampleLomGird(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):
|
||||
|
||||
Reference in New Issue
Block a user