diff --git a/.bumpversion.cfg b/.bumpversion.cfg
index f739d080..ea37cced 100644
--- a/.bumpversion.cfg
+++ b/.bumpversion.cfg
@@ -1,4 +1,4 @@
[bumpversion]
-current_version = 0.1.3
+current_version = 0.1.9
files = setup.py SimPEG/__init__.py docs/conf.py
diff --git a/.travis.yml b/.travis.yml
index c2c672bf..db79b031 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,17 +4,24 @@ python:
sudo: false
+addons:
+ apt:
+ packages:
+ - gcc
+ - gfortran
+ - libopenmpi-dev
+ - libmumps-seq-dev
+ - libblas-dev
+ - liblapack-dev
+
env:
- - TEST_DIR=tests/em/examples
+ - TEST_DIR="tests/mesh tests/base tests/utils"
+ - TEST_DIR=tests/examples
- TEST_DIR=tests/em/fdem/forward
- TEST_DIR=tests/em/fdem/inverse/derivs
- TEST_DIR=tests/em/fdem/inverse/adjoint
- TEST_DIR=tests/em/tdem
- - TEST_DIR=tests/mesh
- TEST_DIR=tests/flow
- - TEST_DIR=tests/utils
- - TEST_DIR=tests/base
- - TEST_DIR=tests/examples
# Setup anaconda
before_install:
@@ -26,9 +33,12 @@ before_install:
# Install packages
install:
- - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython
+ - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose
- pip install nose-cov python-coveralls
- # - pip install -r requirements.txt
+
+ - git clone https://github.com/rowanc1/pymatsolver.git
+ - cd pymatsolver; python setup.py install; cd ..
+
- python setup.py install
- python setup.py build_ext --inplace
diff --git a/SimPEG/EM/Analytics/FDEM.py b/SimPEG/EM/Analytics/FDEM.py
index ce7e623d..9e776fdf 100644
--- a/SimPEG/EM/Analytics/FDEM.py
+++ b/SimPEG/EM/Analytics/FDEM.py
@@ -2,7 +2,6 @@ from __future__ import division
import numpy as np
from scipy.constants import mu_0, pi
from scipy.special import erf
-import matplotlib.pyplot as plt
from SimPEG import Utils
diff --git a/SimPEG/EM/Examples/__init__.py b/SimPEG/EM/Examples/__init__.py
deleted file mode 100644
index eb36678d..00000000
--- a/SimPEG/EM/Examples/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-import CylInversion
diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py
index cce32964..f2167fd8 100644
--- a/SimPEG/EM/FDEM/FDEM.py
+++ b/SimPEG/EM/FDEM/FDEM.py
@@ -66,34 +66,28 @@ class BaseFDEMProblem(BaseEMProblem):
Jv = self.dataPair(self.survey)
for freq in self.survey.freqs:
- dA_du = self.getA(freq) #
- dA_duI = self.Solver(dA_du, **self.solverOpts)
+ A = self.getA(freq) #
+ Ainv = self.Solver(A, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = f[src, ftype]
dA_dm = self.getADeriv_m(freq, u_src, v)
- dRHS_dm = self.getRHSDeriv_m(src, v)
- if dRHS_dm is None:
- du_dm = dA_duI * ( - dA_dm )
- else:
- du_dm = dA_duI * ( - dA_dm + dRHS_dm )
+ dRHS_dm = self.getRHSDeriv_m(freq, src, v)
+ du_dm = Ainv * ( - dA_dm + dRHS_dm )
+
for rx in src.rxList:
- # df_duFun = u.deriv_u(rx.fieldsUsed, m)
df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
- df_du = df_duFun(src, du_dm, adjoint=False)
- if df_du is not None:
- du_dm = df_du
+ df_dudu_dm = df_duFun(src, du_dm, adjoint=False)
df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
df_dm = df_dmFun(src, v, adjoint=False)
- if df_dm is not None:
- du_dm += df_dm
+
+ Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex)
P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m
-
- Jv[src, rx] = P(du_dm)
+ Jv[src, rx] = P(Df_Dm)
return Utils.mkvc(Jv)
@@ -126,30 +120,23 @@ class BaseFDEMProblem(BaseEMProblem):
df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
df_duT = df_duTFun(src, PTv, adjoint=True)
- if df_duT is not None:
- dA_duIT = ATinv * df_duT
- else:
- dA_duIT = ATinv * PTv
+
+ ATinvdf_duT = ATinv * df_duT
- dA_dmT = self.getADeriv_m(freq, u_src, dA_duIT, adjoint=True)
-
- dRHS_dmT = self.getRHSDeriv_m(src, dA_duIT, adjoint=True)
-
- if dRHS_dmT is None:
- du_dmT = - dA_dmT
- else:
- du_dmT = -dA_dmT + dRHS_dmT
+ dA_dmT = self.getADeriv_m(freq, u_src, ATinvdf_duT, adjoint=True)
+ dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True)
+ du_dmT = -dA_dmT + dRHS_dmT
df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
dfT_dm = df_dmFun(src, PTv, adjoint=True)
- if dfT_dm is not None:
- du_dmT += dfT_dm
+
+ du_dmT += dfT_dm
real_or_imag = rx.projComp
- if real_or_imag == 'real':
- Jtv += du_dmT.real
- elif real_or_imag == 'imag':
- Jtv += - du_dmT.real
+ if real_or_imag is 'real':
+ Jtv += np.array(du_dmT,dtype=complex).real
+ elif real_or_imag is 'imag':
+ Jtv += - np.array(du_dmT,dtype=complex).real
else:
raise Exception('Must be real or imag')
@@ -173,10 +160,8 @@ class BaseFDEMProblem(BaseEMProblem):
for i, src in enumerate(Srcs):
smi, sei = src.eval(self)
- if smi is not None:
- S_m[:,i] = Utils.mkvc(smi)
- if sei is not None:
- S_e[:,i] = Utils.mkvc(sei)
+ S_m[:,i] = S_m[:,i] + smi
+ S_e[:,i] = S_e[:,i] + sei
return S_m, S_e
@@ -249,39 +234,21 @@ class Problem_e(BaseFDEMProblem):
C = self.mesh.edgeCurl
MfMui = self.MfMui
- # RHS = C.T * (MfMui * S_m) -1j * omega(freq) * Me * S_e
RHS = C.T * (MfMui * S_m) -1j * omega(freq) * S_e
return RHS
- def getRHSDeriv_m(self, src, v, adjoint=False):
+ def getRHSDeriv_m(self, freq, src, v, adjoint=False):
C = self.mesh.edgeCurl
MfMui = self.MfMui
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
if adjoint:
dRHS = MfMui * (C * v)
- S_mDerivv = S_mDeriv(dRHS)
- S_eDerivv = S_eDeriv(v)
- if S_mDerivv is not None and S_eDerivv is not None:
- return S_mDerivv - 1j * omega(freq) * S_eDerivv
- elif S_mDerivv is not None:
- return S_mDerivv
- elif S_eDerivv is not None:
- return - 1j * omega(freq) * S_eDerivv
- else:
- return None
- else:
- S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v)
+ return S_mDeriv(dRHS) - 1j * omega(freq) * S_eDeriv(v)
- if S_mDerivv is not None and S_eDerivv is not None:
- return C.T * (MfMui * S_mDerivv) -1j * omega(freq) * S_eDerivv
- elif S_mDerivv is not None:
- return C.T * (MfMui * S_mDerivv)
- elif S_eDerivv is not None:
- return -1j * omega(freq) * S_eDerivv
- else:
- return None
+ else:
+ return C.T * (MfMui * S_mDeriv(v)) -1j * omega(freq) * S_eDeriv(v)
class Problem_b(BaseFDEMProblem):
@@ -362,7 +329,6 @@ class Problem_b(BaseFDEMProblem):
S_m, S_e = self.getSourceTerm(freq)
C = self.mesh.edgeCurl
MeSigmaI = self.MeSigmaI
- # Me = self.Me
RHS = S_m + C * ( MeSigmaI * S_e )
@@ -372,51 +338,28 @@ class Problem_b(BaseFDEMProblem):
return RHS
- def getRHSDeriv_m(self, src, v, adjoint=False):
+ def getRHSDeriv_m(self, freq, src, v, adjoint=False):
C = self.mesh.edgeCurl
S_m, S_e = src.eval(self)
MfMui = self.MfMui
- # Me = self.Me
if self._makeASymmetric and adjoint:
v = self.MfMui * v
- if S_e is not None:
- MeSigmaIDeriv = self.MeSigmaIDeriv(S_e)
- if not adjoint:
- RHSderiv = C * (MeSigmaIDeriv * v)
- elif adjoint:
- RHSderiv = MeSigmaIDeriv.T * (C.T * v)
- else:
- RHSderiv = None
-
+ MeSigmaIDeriv = self.MeSigmaIDeriv(S_e)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
- S_mDeriv, S_eDeriv = S_mDeriv(v), S_eDeriv(v)
- if S_mDeriv is not None and S_eDeriv is not None:
- if not adjoint:
- SrcDeriv = S_mDeriv + C * (self.MeSigmaI * S_eDeriv)
- elif adjoint:
- SrcDeriv = S_mDeriv + Self.MeSigmaI.T * ( C.T * S_eDeriv)
- elif S_mDeriv is not None:
- SrcDeriv = S_mDeriv
- elif S_eDeriv is not None:
- if not adjoint:
- SrcDeriv = C * (self.MeSigmaI * S_eDeriv)
- elif adjoint:
- SrcDeriv = self.MeSigmaI.T * ( C.T * S_eDeriv)
- else:
- SrcDeriv = None
- if RHSderiv is not None and SrcDeriv is not None:
- RHSderiv += SrcDeriv
- elif SrcDeriv is not None:
- RHSderiv = SrcDeriv
+ if not adjoint:
+ RHSderiv = C * (MeSigmaIDeriv * v)
+ SrcDeriv = S_mDeriv(v) + C * (self.MeSigmaI * S_eDeriv(v))
+ elif adjoint:
+ RHSderiv = MeSigmaIDeriv.T * (C.T * v)
+ SrcDeriv = S_mDeriv(v) + self.MeSigmaI.T * (C.T * S_eDeriv(v))
- if RHSderiv is not None:
- if self._makeASymmetric is True and not adjoint:
- return MfMui.T * RHSderiv
+ if self._makeASymmetric is True and not adjoint:
+ return MfMui.T * (SrcDeriv + RHSderiv)
- return RHSderiv
+ return RHSderiv + SrcDeriv
@@ -519,7 +462,7 @@ class Problem_j(BaseFDEMProblem):
return RHS
- def getRHSDeriv_m(self, src, v, adjoint=False):
+ def getRHSDeriv_m(self, freq, src, v, adjoint=False):
C = self.mesh.edgeCurl
MeMuI = self.MeMuI
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
@@ -528,27 +471,10 @@ class Problem_j(BaseFDEMProblem):
if self._makeASymmetric:
MfRho = self.MfRho
v = MfRho*v
- S_mDerivv = S_mDeriv(MeMuI.T * (C.T * v))
- S_eDerivv = S_eDeriv(v)
- if S_mDerivv is not None and S_eDerivv is not None:
- return S_mDerivv - 1j * omega(freq) * S_eDerivv
- elif S_mDerivv is not None:
- return S_mDerivv
- elif S_eDerivv is not None:
- return - 1j * omega(freq) * S_eDerivv
- else:
- return None
- else:
- S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v)
+ return S_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * S_eDeriv(v)
- if S_mDerivv is not None and S_eDerivv is not None:
- RHSDeriv = C * (MeMuI * S_mDerivv) - 1j * omega(freq) * S_eDerivv
- elif S_mDerivv is not None:
- RHSDeriv = C * (MeMuI * S_mDerivv)
- elif S_eDerivv is not None:
- RHSDeriv = - 1j * omega(freq) * S_eDerivv
- else:
- return None
+ else:
+ RHSDeriv = C * (MeMuI * S_mDeriv(v)) - 1j * omega(freq) * S_eDeriv(v)
if self._makeASymmetric:
MfRho = self.MfRho
@@ -627,35 +553,18 @@ class Problem_h(BaseFDEMProblem):
return RHS
- def getRHSDeriv_m(self, src, v, adjoint=False):
+ def getRHSDeriv_m(self, freq, src, v, adjoint=False):
_, S_e = src.eval(self)
C = self.mesh.edgeCurl
MfRho = self.MfRho
- RHSDeriv = None
-
- if S_e is not None:
- MfRhoDeriv = self.MfRhoDeriv(S_e)
- if not adjoint:
- RHSDeriv = C.T * (MfRhoDeriv * v)
- elif adjoint:
- RHSDeriv = MfRhoDeriv.T * (C * v)
+ MfRhoDeriv = self.MfRhoDeriv(S_e)
+ if not adjoint:
+ RHSDeriv = C.T * (MfRhoDeriv * v)
+ elif adjoint:
+ RHSDeriv = MfRhoDeriv.T * (C * v)
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
- S_mDeriv = S_mDeriv(v)
- S_eDeriv = S_eDeriv(v)
-
- if S_mDeriv is not None:
- if RHSDeriv is not None:
- RHSDeriv += S_mDeriv(v)
- else:
- RHSDeriv = S_mDeriv(v)
- if S_eDeriv is not None:
- if RHSDeriv is not None:
- RHSDeriv += C.T * (MfRho * S_e)
- else:
- RHSDeriv = C.T * (MfRho * S_e)
-
- return RHSDeriv
+ return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v))
diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py
index d83877c3..8f6fafe9 100644
--- a/SimPEG/EM/FDEM/FieldsFDEM.py
+++ b/SimPEG/EM/FDEM/FieldsFDEM.py
@@ -3,6 +3,7 @@ import scipy.sparse as sp
import SimPEG
from SimPEG import Utils
from SimPEG.EM.Utils import omega
+from SimPEG.Utils import Zero, Identity
class Fields(SimPEG.Problem.Fields):
@@ -32,8 +33,7 @@ class Fields_e(Fields):
ePrimary = np.zeros_like(eSolution)
for i, src in enumerate(srcList):
ep = src.ePrimary(self.prob)
- if ep is not None:
- ePrimary[:,i] = ep
+ ePrimary[:,i] = ePrimary[:,i] + ep
return ePrimary
def _eSecondary(self, eSolution, srcList):
@@ -43,18 +43,17 @@ class Fields_e(Fields):
return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
def _eDeriv_u(self, src, v, adjoint = False):
- return None
+ return Identity()*v
def _eDeriv_m(self, src, v, adjoint = False):
# assuming primary does not depend on the model
- return None
+ return Zero()
def _bPrimary(self, eSolution, srcList):
bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex)
for i, src in enumerate(srcList):
bp = src.bPrimary(self.prob)
- if bp is not None:
- bPrimary[:,i] += bp
+ bPrimary[:,i] = bPrimary[:,i] + bp
return bPrimary
def _bSecondary(self, eSolution, srcList):
@@ -63,8 +62,7 @@ class Fields_e(Fields):
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
S_m, _ = src.eval(self.prob)
- if S_m is not None:
- b[:,i] += 1./(1j*omega(src.freq)) * S_m
+ b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m
return b
def _bSecondaryDeriv_u(self, src, v, adjoint = False):
@@ -76,9 +74,7 @@ class Fields_e(Fields):
def _bSecondaryDeriv_m(self, src, v, adjoint = False):
S_mDeriv, _ = src.evalDeriv(self.prob, adjoint)
S_mDeriv = S_mDeriv(v)
- if S_mDeriv is not None:
- return 1./(1j * omega(src.freq)) * S_mDeriv
- return None
+ return 1./(1j * omega(src.freq)) * S_mDeriv
def _b(self, eSolution, srcList):
return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
@@ -118,8 +114,7 @@ class Fields_b(Fields):
bPrimary = np.zeros_like(bSolution)
for i, src in enumerate(srcList):
bp = src.bPrimary(self.prob)
- if bp is not None:
- bPrimary[:,i] = bp
+ bPrimary[:,i] = bPrimary[:,i] + bp
return bPrimary
def _bSecondary(self, bSolution, srcList):
@@ -129,26 +124,24 @@ class Fields_b(Fields):
return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList)
def _bDeriv_u(self, src, v, adjoint=False):
- return None
+ return Identity()*v
def _bDeriv_m(self, src, v, adjoint=False):
# assuming primary does not depend on the model
- return None
+ return Zero()
def _ePrimary(self, bSolution, srcList):
ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex)
for i,src in enumerate(srcList):
ep = src.ePrimary(self.prob)
- if ep is not None:
- ePrimary[:,i] = ep
+ ePrimary[:,i] = ePrimary[:,i] + ep
return ePrimary
def _eSecondary(self, bSolution, srcList):
e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution))
for i,src in enumerate(srcList):
_,S_e = src.eval(self.prob)
- if S_e is not None:
- e[:,i] += -self._MeSigmaI * S_e
+ e[:,i] = e[:,i]+ -self._MeSigmaI * S_e
return e
def _eSecondaryDeriv_u(self, src, v, adjoint=False):
@@ -166,8 +159,7 @@ class Fields_b(Fields):
Me = Me.T
w = self._edgeCurl.T * (self._MfMui * bSolution)
- if S_e is not None:
- w += -Utils.mkvc(Me * S_e,2)
+ w = w - Utils.mkvc(Me * S_e,2)
if not adjoint:
de_dm = self._MeSigmaIDeriv(w) * v
@@ -177,8 +169,7 @@ class Fields_b(Fields):
_, S_eDeriv = src.evalDeriv(self.prob, adjoint)
Se_Deriv = S_eDeriv(v)
- if Se_Deriv is not None:
- de_dm += -self._MeSigmaI * Se_Deriv
+ de_dm = de_dm - self._MeSigmaI * Se_Deriv
return de_dm
@@ -219,8 +210,7 @@ class Fields_j(Fields):
jPrimary = np.zeros_like(jSolution,dtype = complex)
for i, src in enumerate(srcList):
jp = src.jPrimary(self.prob)
- if jp is not None:
- jPrimary[:,i] += jp
+ jPrimary[:,i] = jPrimary[:,i] + jp
return jPrimary
def _jSecondary(self, jSolution, srcList):
@@ -230,18 +220,17 @@ class Fields_j(Fields):
return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList)
def _jDeriv_u(self, src, v, adjoint=False):
- return None
+ return Identity()*v
def _jDeriv_m(self, src, v, adjoint=False):
# assuming primary does not depend on the model
- return None
+ return Zero()
def _hPrimary(self, jSolution, srcList):
hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex)
for i, src in enumerate(srcList):
hp = src.hPrimary(self.prob)
- if hp is not None:
- hPrimary[:,i] = hp
+ hPrimary[:,i] = hPrimary[:,i] + hp
return hPrimary
def _hSecondary(self, jSolution, srcList):
@@ -249,8 +238,7 @@ class Fields_j(Fields):
for i, src in enumerate(srcList):
h[:,i] *= -1./(1j*omega(src.freq))
S_m,_ = src.eval(self.prob)
- if S_m is not None:
- h[:,i] += 1./(1j*omega(src.freq)) * self._MeMuI * (S_m)
+ h[:,i] = h[:,i]+ 1./(1j*omega(src.freq)) * self._MeMuI * (S_m)
return h
def _hSecondaryDeriv_u(self, src, v, adjoint=False):
@@ -276,12 +264,10 @@ class Fields_j(Fields):
if not adjoint:
S_mDeriv = S_mDeriv(v)
- if S_mDeriv is not None:
- hDeriv_m += 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv)
+ hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv)
elif adjoint:
S_mDeriv = S_mDeriv(Me.T * (MeMuI.T * v))
- if S_mDeriv is not None:
- hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv
+ hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv
return hDeriv_m
@@ -320,9 +306,8 @@ class Fields_h(Fields):
hPrimary = np.zeros_like(hSolution,dtype = complex)
for i, src in enumerate(srcList):
hp = src.hPrimary(self.prob)
- if hp is not None:
- hPrimary[:,i] += hp
- return hPrimary
+ hPrimary[:,i] = hPrimary[:,i] + hp
+ return hPrimary
def _hSecondary(self, hSolution, srcList):
return hSolution
@@ -331,26 +316,24 @@ class Fields_h(Fields):
return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList)
def _hDeriv_u(self, src, v, adjoint=False):
- return None
+ return Identity()*v
def _hDeriv_m(self, src, v, adjoint=False):
# assuming primary does not depend on the model
- return None
+ return Zero()
def _jPrimary(self, hSolution, srcList):
jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]], dtype = complex)
for i, src in enumerate(srcList):
jp = src.jPrimary(self.prob)
- if jp is not None:
- jPrimary[:,i] = jp
+ jPrimary[:,i] = jPrimary[:,i] + jp
return jPrimary
def _jSecondary(self, hSolution, srcList):
j = self._edgeCurl*hSolution
for i, src in enumerate(srcList):
_,S_e = src.eval(self.prob)
- if S_e is not None:
- j[:,i] += -S_e
+ j[:,i] = j[:,i]+ -S_e
return j
def _jSecondaryDeriv_u(self, src, v, adjoint=False):
@@ -362,9 +345,7 @@ class Fields_h(Fields):
def _jSecondaryDeriv_m(self, src, v, adjoint=False):
_,S_eDeriv = src.evalDeriv(self.prob, adjoint)
S_eDeriv = S_eDeriv(v)
- if S_eDeriv is not None:
- return -S_eDeriv
- return None
+ return -S_eDeriv
def _j(self, hSolution, srcList):
return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList)
diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py
index 5dcf9407..b29768ac 100644
--- a/SimPEG/EM/FDEM/SrcFDEM.py
+++ b/SimPEG/EM/FDEM/SrcFDEM.py
@@ -1,6 +1,7 @@
from SimPEG import Survey, Problem, Utils, np, sp
from scipy.constants import mu_0
from SimPEG.EM.Utils import *
+from SimPEG.Utils import Zero
# from SurveyFDEM import Rx
@@ -18,28 +19,28 @@ class BaseSrc(Survey.BaseSrc):
return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
def bPrimary(self, prob):
- return None
+ return Zero()
def hPrimary(self, prob):
- return None
+ return Zero()
def ePrimary(self, prob):
- return None
+ return Zero()
def jPrimary(self, prob):
- return None
+ return Zero()
def S_m(self, prob):
- return None
+ return Zero()
def S_e(self, prob):
- return None
+ return Zero()
def S_mDeriv(self, prob, v, adjoint = False):
- return None
+ return Zero()
def S_eDeriv(self, prob, v, adjoint = False):
- return None
+ return Zero()
class RawVec_e(BaseSrc):
@@ -51,30 +52,14 @@ class RawVec_e(BaseSrc):
:param rxList: receiver list
"""
- def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
+ def __init__(self, rxList, freq, S_e): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
self._S_e = np.array(S_e,dtype=complex)
- self._ePrimary = ePrimary
- self._bPrimary = bPrimary
- self._hPrimary = hPrimary
- self._jPrimary = jPrimary
self.freq = float(freq)
BaseSrc.__init__(self, rxList)
def S_e(self, prob):
return self._S_e
- def ePrimary(self, prob):
- return self._ePrimary
-
- def bPrimary(self, prob):
- return self._bPrimary
-
- def hPrimary(self, prob):
- return self._hPrimary
-
- def jPrimary(self, prob):
- return self._jPrimary
-
class RawVec_m(BaseSrc):
"""
@@ -85,32 +70,16 @@ class RawVec_m(BaseSrc):
:param rxList: receiver list
"""
- def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
+ def __init__(self, rxList, freq, S_m, integrate = True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
self._S_m = np.array(S_m,dtype=complex)
self.freq = float(freq)
self.integrate = integrate
- self._ePrimary = np.array(ePrimary,dtype=complex)
- self._bPrimary = np.array(bPrimary,dtype=complex)
- self._hPrimary = np.array(hPrimary,dtype=complex)
- self._jPrimary = np.array(jPrimary,dtype=complex)
BaseSrc.__init__(self, rxList)
def S_m(self, prob):
return self._S_m
- def ePrimary(self, prob):
- return self._ePrimary
-
- def bPrimary(self, prob):
- return self._bPrimary
-
- def hPrimary(self, prob):
- return self._hPrimary
-
- def jPrimary(self, prob):
- return self._jPrimary
-
class RawVec(BaseSrc):
"""
@@ -192,7 +161,7 @@ class MagDipole(BaseSrc):
def S_e(self, prob):
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
- return None
+ return Zero()
else:
eqLocs = prob._eqLocs
@@ -261,7 +230,7 @@ class MagDipole_Bfield(BaseSrc):
def S_e(self, prob):
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
- return None
+ return Zero()
else:
eqLocs = prob._eqLocs
@@ -329,7 +298,7 @@ class CircularLoop(BaseSrc):
def S_e(self, prob):
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
- return None
+ return Zero()
else:
eqLocs = prob._eqLocs
diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py
index 150a6c00..f60cbfdf 100644
--- a/SimPEG/EM/FDEM/SurveyFDEM.py
+++ b/SimPEG/EM/FDEM/SurveyFDEM.py
@@ -1,8 +1,10 @@
import SimPEG
from SimPEG.EM.Utils import *
from scipy.constants import mu_0
+from SimPEG.Utils import Zero, Identity
import SrcFDEM as Src
+
####################################################
# Receivers
####################################################
diff --git a/SimPEG/EM/Examples/CylInversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py
similarity index 92%
rename from SimPEG/EM/Examples/CylInversion.py
rename to SimPEG/Examples/EM_FDEM_1D_Inversion.py
index cfcfcfc1..aba70f4b 100644
--- a/SimPEG/EM/Examples/CylInversion.py
+++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py
@@ -1,9 +1,16 @@
from SimPEG import *
import SimPEG.EM as EM
from scipy.constants import mu_0
-import matplotlib.pyplot as plt
+
def run(plotIt=True):
+ """
+ EM: FDEM: 1D: Inversion
+ =======================
+
+ Here we will create and run a FDEM 1D inversion.
+
+ """
cs, ncx, ncz, npad = 5., 25, 15, 15
hx = [(cs,ncx), (cs,npad,1.3)]
@@ -24,6 +31,7 @@ def run(plotIt=True):
if plotIt:
+ import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (3, 6))
plt.semilogx(sigma[active], mesh.vectorCCz[active])
ax.set_ylim(-600, 0)
@@ -53,6 +61,7 @@ def run(plotIt=True):
survey.Wd = 1/(abs(survey.dobs)*std)
if plotIt:
+ import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (10, 6))
ax.loglog(rx.times, dtrue, 'b.-')
ax.loglog(rx.times, survey.dobs, 'r.-')
@@ -81,6 +90,7 @@ def run(plotIt=True):
mopt = inv.run(m0)
if plotIt:
+ import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1, figsize = (3, 6))
plt.semilogx(sigma[active], mesh.vectorCCz[active])
plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])
diff --git a/SimPEG/FLOW/Examples/Celia1990.py b/SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
similarity index 53%
rename from SimPEG/FLOW/Examples/Celia1990.py
rename to SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
index 24ae82a6..4c90ea7b 100644
--- a/SimPEG/FLOW/Examples/Celia1990.py
+++ b/SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
@@ -1,8 +1,40 @@
from SimPEG import *
from SimPEG.FLOW import Richards
-import matplotlib.pyplot as plt
def run(plotIt=True):
+ """
+ FLOW: Richards: 1D: Celia1990
+ =============================
+
+ There are two different forms of Richards equation that differ
+ on how they deal with the non-linearity in the time-stepping term.
+
+ The most fundamental form, referred to as the
+ 'mixed'-form of Richards Equation Celia1990_
+
+ .. math::
+
+ \\frac{\partial \\theta(\psi)}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0
+ \quad \psi \in \Omega
+
+ where \\\\(\\\\theta\\\\) is water content, and \\\\(\\\\psi\\\\) is pressure head.
+ This formulation of Richards equation is called the
+ 'mixed'-form because the equation is parameterized in \\\\(\\\\psi\\\\)
+ but the time-stepping is in terms of \\\\(\\\\theta\\\\).
+
+ As noted in Celia1990_ the 'head'-based form of Richards
+ equation can be written in the continuous form as:
+
+ .. math::
+
+ \\frac{\partial \\theta}{\partial \psi}\\frac{\partial \psi}{\partial t} - \\nabla \cdot k(\psi) \\nabla \psi - \\frac{\partial k(\psi)}{\partial z} = 0 \quad \psi \in \Omega
+
+ However, it can be shown that this does not conserve mass in the discrete formulation.
+
+ Here we reproduce the results from Celia1990_ demonstrating the head-based formulation and the mixed-formulation.
+
+ .. _Celia1990: http://www.webpages.uidaho.edu/ch/papers/Celia.pdf
+ """
M = Mesh.TensorMesh([np.ones(40)])
M.setCellGradBC('dirichlet')
params = Richards.Empirical.HaverkampParams().celia1990
@@ -28,6 +60,7 @@ def run(plotIt=True):
Hs_H120= getFields(120.,'head')
if not plotIt:return
+ import matplotlib.pyplot as plt
plt.figure(figsize=(13,5))
plt.subplot(121)
plt.plot(40-M.gridCC, Hs_M10[-1],'b-')
@@ -47,6 +80,7 @@ def run(plotIt=True):
plt.xlabel('Depth, cm')
plt.ylabel('Pressure Head, cm')
plt.legend(('$\Delta t$ = 10 sec','$\Delta t$ = 30 sec','$\Delta t$ = 120 sec'))
+ plt.show()
if __name__ == '__main__':
run()
diff --git a/SimPEG/Examples/DCfwd.py b/SimPEG/Examples/Forward_BasicDirectCurrent.py
similarity index 95%
rename from SimPEG/Examples/DCfwd.py
rename to SimPEG/Examples/Forward_BasicDirectCurrent.py
index 33e7aad5..efbb287c 100644
--- a/SimPEG/Examples/DCfwd.py
+++ b/SimPEG/Examples/Forward_BasicDirectCurrent.py
@@ -1,7 +1,4 @@
from SimPEG import Mesh, Utils, np, SolverLU
-import matplotlib.pyplot as plt
-import matplotlib
-from matplotlib.mlab import griddata
## 2D DC forward modeling example with Tensor and Curvilinear Meshes
@@ -12,7 +9,6 @@ def run(plotIt=True):
tM = Mesh.TensorMesh(sz)
# Curvilinear Mesh
rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
-
# Step2: Direct Current (DC) operator
def DCfun(mesh, pts):
D = mesh.faceDiv
@@ -39,6 +35,11 @@ def run(plotIt=True):
phirM = AinvrM*rhsrM
if not plotIt: return
+
+ import matplotlib.pyplot as plt
+ import matplotlib
+ from matplotlib.mlab import griddata
+
#Step4: Making Figure
fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
label = ["(a)", "(b)"]
@@ -69,6 +70,7 @@ def run(plotIt=True):
else:
axes[i].set_ylabel(" ")
axes[i].set_xlabel("x")
+ plt.show()
if __name__ == '__main__':
diff --git a/SimPEG/Examples/Linear.py b/SimPEG/Examples/Inversion_Linear.py
similarity index 50%
rename from SimPEG/Examples/Linear.py
rename to SimPEG/Examples/Inversion_Linear.py
index b065682b..a8a0eddc 100644
--- a/SimPEG/Examples/Linear.py
+++ b/SimPEG/Examples/Inversion_Linear.py
@@ -1,29 +1,39 @@
from SimPEG import *
-class LinearSurvey(Survey.BaseSurvey):
- def projectFields(self, u):
- return u
-class LinearProblem(Problem.BaseProblem):
- """docstring for LinearProblem"""
+def run(N=100, plotIt=True):
+ """
+ Inversion: Linear Problem
+ =========================
- surveyPair = LinearSurvey
+ Here we go over the basics of creating a linear problem and inversion.
- def __init__(self, mesh, G, **kwargs):
- Problem.BaseProblem.__init__(self, mesh, **kwargs)
- self.G = G
+ """
- def fields(self, m, u=None):
- return self.G.dot(m)
+ class LinearSurvey(Survey.BaseSurvey):
+ def projectFields(self, u):
+ return u
- def Jvec(self, m, v, u=None):
- return self.G.dot(v)
+ class LinearProblem(Problem.BaseProblem):
- def Jtvec(self, m, v, u=None):
- return self.G.T.dot(v)
+ surveyPair = LinearSurvey
+
+ def __init__(self, mesh, G, **kwargs):
+ Problem.BaseProblem.__init__(self, mesh, **kwargs)
+ self.G = G
+
+ def fields(self, m, u=None):
+ return self.G.dot(m)
+
+ def Jvec(self, m, v, u=None):
+ return self.G.dot(v)
+
+ def Jtvec(self, m, v, u=None):
+ return self.G.T.dot(v)
-def run(N, plotIt=True):
+ np.random.seed(1)
+
mesh = Mesh.TensorMesh([N])
nk = 20
@@ -52,7 +62,7 @@ def run(N, plotIt=True):
reg = Regularization.Tikhonov(mesh)
dmis = DataMisfit.l2_DataMisfit(survey)
- opt = Optimization.InexactGaussNewton(maxIter=20)
+ opt = Optimization.InexactGaussNewton(maxIter=35)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
beta = Directives.BetaSchedule()
betaest = Directives.BetaEstimate_ByEig()
@@ -63,16 +73,18 @@ def run(N, plotIt=True):
if plotIt:
import matplotlib.pyplot as plt
- plt.figure(1)
- for i in range(prob.G.shape[0]):
- plt.plot(prob.G[i,:])
- plt.figure(2)
- plt.plot(M.vectorCCx, survey.mtrue, 'b-')
- plt.plot(M.vectorCCx, mrec, 'r-')
+ fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
+ for i in range(prob.G.shape[0]):
+ axes[0].plot(prob.G[i,:])
+ axes[0].set_title('Columns of matrix G')
+
+ axes[1].plot(M.vectorCCx, survey.mtrue, 'b-')
+ axes[1].plot(M.vectorCCx, mrec, 'r-')
+ axes[1].legend(('True Model', 'Recovered Model'))
plt.show()
return prob, survey, mesh, mrec
if __name__ == '__main__':
- run(100)
+ run()
diff --git a/SimPEG/Examples/Mesh_Basic_PlotImage.py b/SimPEG/Examples/Mesh_Basic_PlotImage.py
new file mode 100644
index 00000000..3154f281
--- /dev/null
+++ b/SimPEG/Examples/Mesh_Basic_PlotImage.py
@@ -0,0 +1,46 @@
+from SimPEG import *
+
+def run(plotIt=True):
+ """
+ Mesh: Basic: PlotImage
+ ======================
+
+ You can use M.PlotImage to plot images on all of the Meshes.
+
+
+ """
+ M = Mesh.TensorMesh([32,32])
+ v = Utils.ModelBuilder.randomModel(M.vnC, seed=789)
+ v = Utils.mkvc(v)
+
+ O = Mesh.TreeMesh([32,32])
+ O.refine(1)
+ def function(cell):
+ if (cell.center[0] < 0.75 and cell.center[0] > 0.25 and
+ cell.center[1] < 0.75 and cell.center[1] > 0.25):return 5
+ if (cell.center[0] < 0.9 and cell.center[0] > 0.1 and
+ cell.center[1] < 0.9 and cell.center[1] > 0.1):return 4
+ return 3
+ O.refine(function)
+
+ P = M.getInterpolationMat(O.gridCC, 'CC')
+
+ ov = P * v
+
+ if plotIt:
+ import matplotlib.pyplot as plt
+
+ fig, axes = plt.subplots(1,2,figsize=(10,5))
+
+ out = M.plotImage(v, grid=True, ax=axes[0])
+ cb = plt.colorbar(out[0], ax=axes[0]); cb.set_label("Random Field")
+ axes[0].set_title('TensorMesh')
+
+ out = O.plotImage(ov, grid=True, ax=axes[1], clim=[0,1])
+ cb = plt.colorbar(out[0], ax=axes[1]); cb.set_label("Random Field")
+ axes[1].set_title('TreeMesh')
+
+ plt.show()
+
+if __name__ == '__main__':
+ run()
diff --git a/SimPEG/Examples/Mesh_Basic_Types.py b/SimPEG/Examples/Mesh_Basic_Types.py
new file mode 100644
index 00000000..430fe698
--- /dev/null
+++ b/SimPEG/Examples/Mesh_Basic_Types.py
@@ -0,0 +1,30 @@
+from SimPEG import *
+
+def run(plotIt=True):
+ """
+ Mesh: Basic: Types
+ ==================
+
+ Here we show SimPEG used to create three different types of meshes.
+
+ """
+ sz = [16,16]
+ tM = Mesh.TensorMesh(sz)
+ qM = Mesh.TreeMesh(sz)
+ qM.refine(lambda cell: 4 if np.sqrt(((np.r_[cell.center]-0.5)**2).sum()) < 0.4 else 3)
+ rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
+
+ if plotIt:
+ import matplotlib.pyplot as plt
+ fig, axes = plt.subplots(1,3,figsize=(14,5))
+ opts = {}
+ tM.plotGrid(ax=axes[0], **opts)
+ axes[0].set_title('TensorMesh')
+ qM.plotGrid(ax=axes[1], **opts)
+ axes[1].set_title('TreeMesh')
+ rM.plotGrid(ax=axes[2], **opts)
+ axes[2].set_title('CurvilinearMesh')
+ plt.show()
+
+if __name__ == '__main__':
+ run()
diff --git a/SimPEG/Examples/Mesh_Operators_CahnHilliard.py b/SimPEG/Examples/Mesh_Operators_CahnHilliard.py
new file mode 100644
index 00000000..8e42e617
--- /dev/null
+++ b/SimPEG/Examples/Mesh_Operators_CahnHilliard.py
@@ -0,0 +1,105 @@
+from SimPEG import *
+
+def run(plotIt=True, n=60):
+ """
+ Mesh: Operators: Cahn Hilliard
+ ==============================
+
+ This example is based on the example in the FiPy_ library.
+ Please see their documentation for more information about the Cahn-Hilliard equation.
+
+ The "Cahn-Hilliard" equation separates a field \\\\( \\\\phi \\\\) into 0 and 1 with smooth transitions.
+
+ .. math::
+
+ \\frac{\partial \phi}{\partial t} = \\nabla \cdot D \\nabla \left( \\frac{\partial f}{\partial \phi} - \epsilon^2 \\nabla^2 \phi \\right)
+
+ Where \\\\( f \\\\) is the energy function \\\\( f = ( a^2 / 2 )\\\\phi^2(1 - \\\\phi)^2 \\\\)
+ which drives \\\\( \\\\phi \\\\) towards either 0 or 1, this competes with the term
+ \\\\(\\\\epsilon^2 \\\\nabla^2 \\\\phi \\\\) which is a diffusion term that creates smooth changes in \\\\( \\\\phi \\\\).
+ The equation can be factored:
+
+ .. math::
+
+ \\frac{\partial \phi}{\partial t} = \\nabla \cdot D \\nabla \psi \\\\
+ \psi = \\frac{\partial^2 f}{\partial \phi^2} (\phi - \phi^{\\text{old}}) + \\frac{\partial f}{\partial \phi} - \epsilon^2 \\nabla^2 \phi
+
+ Here we will need the derivatives of \\\\( f \\\\):
+
+ .. math::
+
+ \\frac{\partial f}{\partial \phi} = (a^2/2)2\phi(1-\phi)(1-2\phi)
+ \\frac{\partial^2 f}{\partial \phi^2} = (a^2/2)2[1-6\phi(1-\phi)]
+
+ The implementation below uses backwards Euler in time with an exponentially increasing time step.
+ The initial \\\\( \\\\phi \\\\) is a normally distributed field with a standard deviation of 0.1 and mean of 0.5.
+ The grid is 60x60 and takes a few seconds to solve ~130 times. The results are seen below, and you can see the
+ field separating as the time increases.
+
+ .. _FiPy: http://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2DCoupled.html
+
+ """
+
+ np.random.seed(5)
+
+ # Here we are going to rearrange the equations:
+
+ # (phi_ - phi)/dt = A*(d2fdphi2*(phi_ - phi) + dfdphi - L*phi_)
+ # (phi_ - phi)/dt = A*(d2fdphi2*phi_ - d2fdphi2*phi + dfdphi - L*phi_)
+ # (phi_ - phi)/dt = A*d2fdphi2*phi_ + A*( - d2fdphi2*phi + dfdphi - L*phi_)
+ # phi_ - phi = dt*A*d2fdphi2*phi_ + dt*A*(- d2fdphi2*phi + dfdphi - L*phi_)
+ # phi_ - dt*A*d2fdphi2 * phi_ = dt*A*(- d2fdphi2*phi + dfdphi - L*phi_) + phi
+ # (I - dt*A*d2fdphi2) * phi_ = dt*A*(- d2fdphi2*phi + dfdphi - L*phi_) + phi
+ # (I - dt*A*d2fdphi2) * phi_ = dt*A*dfdphi - dt*A*d2fdphi2*phi - dt*A*L*phi_ + phi
+ # (dt*A*d2fdphi2 - I) * phi_ = dt*A*d2fdphi2*phi + dt*A*L*phi_ - phi - dt*A*dfdphi
+ # (dt*A*d2fdphi2 - I - dt*A*L) * phi_ = (dt*A*d2fdphi2 - I)*phi - dt*A*dfdphi
+
+ h = [(0.25,n)]
+ M = Mesh.TensorMesh([h,h])
+
+ # Constants
+ D = a = epsilon = 1.
+ I = Utils.speye(M.nC)
+
+ # Operators
+ A = D * M.faceDiv * M.cellGrad
+ L = epsilon**2 * M.faceDiv * M.cellGrad
+
+ duration = 75
+ elapsed = 0.
+ dexp = -5
+ phi = np.random.normal(loc=0.5,scale=0.01,size=M.nC)
+ ii, jj = 0, 0
+ PHIS = []
+ capture = np.logspace(-1,np.log10(duration),8)
+ while elapsed < duration:
+ dt = min(100, np.exp(dexp))
+ elapsed += dt
+ dexp += 0.05
+
+ dfdphi = a**2 * 2 * phi * (1 - phi) * (1 - 2 * phi)
+ d2fdphi2 = Utils.sdiag(a**2 * 2 * (1 - 6 * phi * (1 - phi)))
+
+ MAT = (dt*A*d2fdphi2 - I - dt*A*L)
+ rhs = (dt*A*d2fdphi2 - I)*phi - dt*A*dfdphi
+ phi = Solver(MAT)*rhs
+
+ if elapsed > capture[jj]:
+ PHIS += [(elapsed, phi.copy())]
+ jj += 1
+ if ii % 10 == 0: print ii, elapsed
+ ii += 1
+
+ if plotIt:
+ import matplotlib.pyplot as plt
+ fig, axes = plt.subplots(2,4,figsize=(14,6))
+ axes = np.array(axes).flatten().tolist()
+ for ii, ax in zip(np.linspace(0,len(PHIS)-1,len(axes)),axes):
+ ii = int(ii)
+ out = M.plotImage(PHIS[ii][1],ax=ax)
+ ax.axis('off')
+ ax.set_title('Elapsed Time: %4.1f'%PHIS[ii][0])
+ plt.show()
+
+if __name__ == '__main__':
+ run()
diff --git a/SimPEG/Examples/Mesh_QuadTree_Creation.py b/SimPEG/Examples/Mesh_QuadTree_Creation.py
new file mode 100644
index 00000000..ede69a63
--- /dev/null
+++ b/SimPEG/Examples/Mesh_QuadTree_Creation.py
@@ -0,0 +1,28 @@
+from SimPEG import *
+
+def run(plotIt=True):
+ """
+ Mesh: QuadTree: Creation
+ ========================
+
+ You can give the refine method a function, which is evaluated on every cell
+ of the TreeMesh.
+
+ Occasionally it is useful to initially refine to a constant level
+ (e.g. 3 in this 32x32 mesh). This means the function is first evaluated
+ on an 8x8 mesh (2^3).
+
+ """
+ M = Mesh.TreeMesh([32,32])
+ M.refine(3)
+ def function(cell):
+ xyz = cell.center
+ for i in range(3):
+ if np.abs(np.sin(xyz[0]*np.pi*2)*0.5 + 0.5 - xyz[1]) < 0.2*i:
+ return 6-i
+ return 0
+ M.refine(function);
+ if plotIt: M.plotGrid(showIt=True)
+
+if __name__ == '__main__':
+ run()
diff --git a/SimPEG/Examples/Mesh_QuadTree_FaceDiv.py b/SimPEG/Examples/Mesh_QuadTree_FaceDiv.py
new file mode 100644
index 00000000..5bd67929
--- /dev/null
+++ b/SimPEG/Examples/Mesh_QuadTree_FaceDiv.py
@@ -0,0 +1,49 @@
+from SimPEG import *
+
+def run(plotIt=True, n=60):
+ """
+ Mesh: QuadTree: FaceDiv
+ =======================
+
+
+
+ """
+
+
+ M = Mesh.TreeMesh([[(1,16)],[(1,16)]], levels=4)
+ M._refineCell([0,0,0])
+ M._refineCell([0,0,1])
+ M._refineCell([4,4,2])
+ M.__dirty__ = True
+ M.number()
+
+
+ if plotIt:
+ import matplotlib.pyplot as plt
+ fig, axes = plt.subplots(2,1,figsize=(10,10))
+
+ M.plotGrid(cells=True, nodes=False, ax=axes[0])
+ axes[0].axis('off')
+ axes[0].set_title('Simple QuadTree Mesh')
+ axes[0].set_xlim([-1,17])
+ axes[0].set_ylim([-1,17])
+
+ for ii, loc in zip(range(M.nC),M.gridCC):
+ axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='r')
+
+ axes[0].plot(M.gridFx[:,0],M.gridFx[:,1], 'g>')
+ for ii, loc in zip(range(M.nFx),M.gridFx):
+ axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='g')
+
+ axes[0].plot(M.gridFy[:,0],M.gridFy[:,1], 'm^')
+ for ii, loc in zip(range(M.nFy),M.gridFy):
+ axes[0].text(loc[0]+0.2,loc[1]+0.2,'%d'%(ii+M.nFx), color='m')
+
+ axes[1].spy(M.faceDiv)
+ axes[1].set_title('Face Divergence')
+ axes[1].set_ylabel('Cell Number')
+ axes[1].set_xlabel('Face Number')
+ plt.show()
+
+if __name__ == '__main__':
+ run()
diff --git a/SimPEG/Examples/Mesh_QuadTree_HangingNodes.py b/SimPEG/Examples/Mesh_QuadTree_HangingNodes.py
new file mode 100644
index 00000000..11726995
--- /dev/null
+++ b/SimPEG/Examples/Mesh_QuadTree_HangingNodes.py
@@ -0,0 +1,32 @@
+from SimPEG import *
+
+def run(plotIt=True):
+ """
+ Mesh: QuadTree: Hanging Nodes
+ =============================
+
+ You can give the refine method a function, which is evaluated on every cell
+ of the TreeMesh.
+
+ Occasionally it is useful to initially refine to a constant level
+ (e.g. 3 in this 32x32 mesh). This means the function is first evaluated
+ on an 8x8 mesh (2^3).
+
+ """
+ M = Mesh.TreeMesh([8,8])
+ def function(cell):
+ xyz = cell.center
+ dist = ((xyz - [0.25,0.25])**2).sum()**0.5
+ if dist < 0.25:
+ return 3
+ return 2
+ M.refine(function);
+ M.number()
+ if plotIt:
+ import matplotlib.pyplot as plt
+ M.plotGrid(nodes=True, cells=True, facesX=True)
+ plt.legend(('Grid', 'Cell Centers', 'Nodes', 'Hanging Nodes', 'X faces', 'Hanging X faces'))
+ plt.show()
+
+if __name__ == '__main__':
+ run()
diff --git a/SimPEG/Examples/Mesh_Tensor_Creation.py b/SimPEG/Examples/Mesh_Tensor_Creation.py
new file mode 100644
index 00000000..31ad3d69
--- /dev/null
+++ b/SimPEG/Examples/Mesh_Tensor_Creation.py
@@ -0,0 +1,35 @@
+from SimPEG import *
+
+def run(plotIt=True):
+ """
+
+ Mesh: Tensor: Creation
+ ======================
+
+ For tensor meshes, there are some functions that can come
+ in handy. For example, creating mesh tensors can be a bit time
+ consuming, these can be created speedily by just giving numbers
+ and sizes of padding. See the example below, that follows this
+ notation::
+
+ h1 = (
+ (cellSize, numPad, [, increaseFactor]),
+ (cellSize, numCore),
+ (cellSize, numPad, [, increaseFactor])
+ )
+
+ .. note::
+
+ You can center your mesh by passing a 'C' for the x0[i] position.
+ A 'N' will make the entire mesh negative, and a '0' (or a 0) will
+ make the mesh start at zero.
+
+ """
+ h1 = [(10, 5, -1.3), (5, 20), (10, 3, 1.3)]
+ M = Mesh.TensorMesh([h1, h1], x0='CN')
+ if plotIt:
+ M.plotGrid(showIt=True)
+
+if __name__ == '__main__':
+ run()
+
diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py
index 3fec0b99..8431e4ba 100644
--- a/SimPEG/Examples/__init__.py
+++ b/SimPEG/Examples/__init__.py
@@ -1 +1,103 @@
-import Linear, DCfwd
+# Run this file to add imports.
+
+##### AUTOIMPORTS #####
+import EM_FDEM_1D_Inversion
+import FLOW_Richards_1D_Celia1990
+import Forward_BasicDirectCurrent
+import Inversion_Linear
+import Mesh_Basic_PlotImage
+import Mesh_Basic_Types
+import Mesh_Operators_CahnHilliard
+import Mesh_QuadTree_Creation
+import Mesh_QuadTree_FaceDiv
+import Mesh_QuadTree_HangingNodes
+import Mesh_Tensor_Creation
+
+__examples__ = ["EM_FDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"]
+
+##### AUTOIMPORTS #####
+
+if __name__ == '__main__':
+ """
+
+ Run the following to create the examples documentation and add to the imports at the top.
+
+ """
+
+ import shutil, os
+ from SimPEG import Examples
+
+ # Create the examples dir in the docs folder.
+ docExamplesDir = os.path.sep.join(os.path.realpath(__file__).split(os.path.sep)[:-3] + ['docs', 'examples'])
+ shutil.rmtree(docExamplesDir)
+ os.makedirs(docExamplesDir)
+
+ # Get all the python examples in this folder
+ thispath = os.path.sep.join(__file__.split(os.path.sep)[:-1])
+ exfiles = [f[:-3] for f in os.listdir(thispath) if os.path.isfile(os.path.join(thispath, f)) and f.endswith('.py') and not f.startswith('_')]
+
+ # Add the imports to the top in the AUTOIMPORTS section
+ f = file(__file__, 'r')
+ inimports = False
+ out = ''
+ for line in f:
+ if not inimports:
+ out += line
+
+ if line == "##### AUTOIMPORTS #####\n":
+ inimports = not inimports
+ if inimports:
+ out += '\n'.join(["import %s"%_ for _ in exfiles])
+ out += '\n\n__examples__ = ["' + '", "'.join(exfiles)+ '"]\n'
+ out += '\n##### AUTOIMPORTS #####\n'
+ f.close()
+
+ f = file(__file__, 'w')
+ f.write(out)
+ f.close()
+
+
+ def _makeExample(filePath, runFunction):
+ """Makes the example given a path of the file and the run function."""
+ filePath = os.path.realpath(filePath)
+ name = filePath.split(os.path.sep)[-1].rstrip('.pyc').rstrip('.py')
+
+ docstr = runFunction.__doc__
+ if docstr is None:
+ doc = '%s\n%s'%(name.replace('_',' '),'='*len(name))
+ else:
+ doc = '\n'.join([_[8:].rstrip() for _ in docstr.split('\n')])
+
+ out = """.. _examples_%s:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+%s
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.%s.run()
+
+.. literalinclude:: ../../SimPEG/Examples/%s.py
+ :language: python
+ :linenos:
+"""%(name,doc,name,name)
+
+ rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst']))
+
+ print 'Creating: %s.rst'%name
+ f = open(rst, 'w')
+ f.write(out)
+ f.close()
+
+ for ex in dir(Examples):
+ if ex.startswith('_'): continue
+ E = getattr(Examples,ex)
+ _makeExample(E.__file__, E.run)
diff --git a/SimPEG/FLOW/Examples/__init__.py b/SimPEG/FLOW/Examples/__init__.py
deleted file mode 100644
index 7a894e11..00000000
--- a/SimPEG/FLOW/Examples/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-import Celia1990
diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py
index 1e6c91c0..1d5c8c8c 100644
--- a/SimPEG/Mesh/TreeMesh.py
+++ b/SimPEG/Mesh/TreeMesh.py
@@ -90,13 +90,14 @@
#
from SimPEG import np, sp, Utils, Solver
-import matplotlib.pyplot as plt
-import matplotlib
-from mpl_toolkits.mplot3d import Axes3D
-import matplotlib.colors as colors
-import matplotlib.cm as cmx
-import TreeUtils
+try:
+ import TreeUtils
+ _IMPORT_TREEUTILS = True
+except Exception, e:
+ _IMPORT_TREEUTILS = False
+
+
from InnerProducts import InnerProducts
from TensorMesh import TensorMesh, BaseTensorMesh
import time
@@ -108,6 +109,8 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
_meshType = 'TREE'
def __init__(self, h, x0=None, levels=None):
+ if not _IMPORT_TREEUTILS:
+ raise Exception('Could not import the Cython code to run the TreeMesh Try:.\n\npython setup.py build_ext --inplace')
assert type(h) is list, 'h must be a list'
assert len(h) in [2,3], "There is only support for TreeMesh in 2D or 3D."
@@ -1960,11 +1963,18 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
def plotGrid(self, ax=None, showIt=False,
grid=True,
- cells=True, cellLine=False,
+ cells=False, cellLine=False,
nodes=False,
facesX=False, facesY=False, facesZ=False,
edgesX=False, edgesY=False, edgesZ=False):
+
+ import matplotlib.pyplot as plt
+ import matplotlib
+ from mpl_toolkits.mplot3d import Axes3D
+ import matplotlib.colors as colors
+ import matplotlib.cm as cmx
+
# self.number()
axOpts = {'projection':'3d'} if self.dim == 3 else {}
@@ -1975,24 +1985,28 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
fig = ax.figure
if grid:
+ X, Y, Z = [], [], []
for ind in self._sortedCells:
p = self._asPointer(ind)
n = self._cellN(p)
h = self._cellH(p)
- x = [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0]]
- y = [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1]]
if self.dim == 2:
- ax.plot(x,y, 'b-')
+ X += [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0], np.nan]
+ Y += [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1], np.nan]
elif self.dim == 3:
- ax.plot(x,y, 'b-', zs=[n[2]]*5)
- z = [n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2]]
- ax.plot(x,y, 'b-', zs=z)
+ X += [n[0] , n[0] + h[0], n[0] + h[0], n[0] , n[0], np.nan]*2
+ Y += [n[1] , n[1] , n[1] + h[1], n[1] + h[1], n[1], np.nan]*2
+ Z += [n[2]]*5+[np.nan]
+ Z += [n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], n[2] + h[2], np.nan]
sides = [0,0], [h[0],0], [0,h[1]], [h[0],h[1]]
for s in sides:
- x = [n[0] + s[0], n[0] + s[0]]
- y = [n[1] + s[1], n[1] + s[1]]
- z = [n[2] , n[2] + h[2]]
- ax.plot(x,y, 'b-', zs=z)
+ X += [n[0] + s[0], n[0] + s[0]]
+ Y += [n[1] + s[1], n[1] + s[1]]
+ Z += [n[2] , n[2] + h[2]]
+ if self.dim == 2:
+ ax.plot(X,Y, 'b-')
+ elif self.dim == 3:
+ ax.plot(X,Y, 'b-', zs=Z)
if self.dim == 2:
if cells:
@@ -2004,11 +2018,13 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
ax.plot(self._gridN[:,0], self._gridN[:,1], 'ms')
ax.plot(self._gridN[self._hangingN.keys(),0], self._gridN[self._hangingN.keys(),1], 'ms', ms=10, mfc='none', mec='m')
if facesX:
- ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g')
ax.plot(self._gridFx[:,0], self._gridFx[:,1], 'g>')
+ ax.plot(self._gridFx[self._hangingFx.keys(),0], self._gridFx[self._hangingFx.keys(),1], 'gs', ms=10, mfc='none', mec='g')
if facesY:
- ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g')
ax.plot(self._gridFy[:,0], self._gridFy[:,1], 'g^')
+ ax.plot(self._gridFy[self._hangingFy.keys(),0], self._gridFy[self._hangingFy.keys(),1], 'gs', ms=10, mfc='none', mec='g')
+ ax.set_xlabel('x1')
+ ax.set_ylabel('x2')
elif self.dim == 3:
if cells:
ax.plot(self.gridCC[:,0], self.gridCC[:,1], 'r.', zs=self.gridCC[:,2])
@@ -2056,7 +2072,6 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
ind = [key, hf[0]]
ax.plot(self._gridEx[ind,0], self._gridEx[ind,1], 'k:', zs=self._gridEx[ind,2])
-
if edgesY:
ax.plot(self._gridEy[:,0], self._gridEy[:,1], 'k<', zs=self._gridEy[:,2])
ax.plot(self._gridEy[self._hangingEy.keys(),0], self._gridEy[self._hangingEy.keys(),1], 'ks', ms=10, mfc='none', mec='k', zs=self._gridEy[self._hangingEy.keys(),2])
@@ -2072,15 +2087,28 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
for hf in self._hangingEz[key]:
ind = [key, hf[0]]
ax.plot(self._gridEz[ind,0], self._gridEz[ind,1], 'k:', zs=self._gridEz[ind,2])
-
+ ax.set_xlabel('x1')
+ ax.set_ylabel('x2')
+ ax.set_zlabel('x3')
+ ax.grid(True)
if showIt:plt.show()
- def plotImage(self, I, ax=None, showIt=True, grid=False):
+ def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
if self.dim == 3: raise Exception('Use plot slice?')
+
+ import matplotlib.pyplot as plt
+ import matplotlib
+ from mpl_toolkits.mplot3d import Axes3D
+ import matplotlib.colors as colors
+ import matplotlib.cm as cmx
+
if ax is None: ax = plt.subplot(111)
jet = cm = plt.get_cmap('jet')
- cNorm = colors.Normalize(vmin=I.min(), vmax=I.max())
+ cNorm = colors.Normalize(
+ vmin=I.min() if clim is None else clim[0],
+ vmax=I.max() if clim is None else clim[1])
+
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
ax.set_xlim((self.x0[0], self.h[0].sum()))
ax.set_ylim((self.x0[1], self.h[1].sum()))
@@ -2089,8 +2117,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
ax.add_patch(plt.Rectangle((x0[0], x0[1]), sz[0], sz[1], facecolor=scalarMap.to_rgba(I[ii]), edgecolor='k' if grid else 'none'))
# if text: ax.text(self.center[0],self.center[1],self.num)
scalarMap._A = [] # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
- plt.colorbar(scalarMap)
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
if showIt: plt.show()
+ return [scalarMap]
def plotSlice(self, v, vType='CC',
normal='Z', ind=None, grid=True, view='real',
@@ -2102,6 +2132,13 @@ class TreeMesh(BaseTensorMesh, InnerProducts):
assert vType in ['CC','F','E']
assert self.dim == 3
+
+ import matplotlib.pyplot as plt
+ import matplotlib
+ from mpl_toolkits.mplot3d import Axes3D
+ import matplotlib.colors as colors
+ import matplotlib.cm as cmx
+
szSliceDim = len(getattr(self, 'h'+normal.lower())) #: Size of the sliced dimension
if ind is None: ind = int(szSliceDim/2)
assert type(ind) in [int, long], 'ind must be an integer'
@@ -2191,7 +2228,7 @@ class Cell(object):
@property
def center(self):
if getattr(self, '_center', None) is None:
- self._center = self.mesh._cellC(self._pointer)
+ self._center = np.array(self.mesh._cellC(self._pointer))
return self._center
@property
def h(self): return self.mesh._cellH(self._pointer)
@@ -2248,6 +2285,13 @@ class CellLookUpException(TreeException):
if __name__ == '__main__':
+
+ import matplotlib.pyplot as plt
+ import matplotlib
+ from mpl_toolkits.mplot3d import Axes3D
+ import matplotlib.colors as colors
+ import matplotlib.cm as cmx
+
def topo(x):
return np.sin(x*(2.*np.pi))*0.3 + 0.5
diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py
index b078eb66..089d7d9a 100644
--- a/SimPEG/Mesh/View.py
+++ b/SimPEG/Mesh/View.py
@@ -1,8 +1,11 @@
import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib
-from mpl_toolkits.mplot3d import Axes3D
-from SimPEG.Utils import mkvc, animate
+from SimPEG.Utils import mkvc
+try:
+ import matplotlib.pyplot as plt
+ import matplotlib
+ from mpl_toolkits.mplot3d import Axes3D
+except ImportError, e:
+ print 'Trouble importing matplotlib.'
class TensorView(object):
@@ -479,63 +482,6 @@ class TensorView(object):
ax.grid(True)
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'
- if ax is None: ax = plt.subplot(111)
- I = mesh.r(var,'CC','CC','M')
- axes = [p for p in 'xyz' if p not in normal.lower()]
- if normal is 'x': I = I[index,:,:]
- if normal is 'y': I = I[:,index,:]
- if normal is 'z': I = I[:,:,index]
- if clim is None: clim = [I.min(),I.max()]
- p = ax.pcolormesh(getattr(mesh,'vectorN'+axes[0]),getattr(mesh,'vectorN'+axes[1]),I.T,vmin=clim[0],vmax=clim[1])
- ax.axis('tight')
- ax.set_xlabel(axes[0])
- ax.set_ylabel(axes[1])
- return p
-
- def videoSlicer(mesh,var,imageType='CC',normal='z',figsize=(10,8)):
- assert mesh.dim > 2, 'This is for 3D meshes only.'
- # First set up the figure, the axis, and the plot element we want to animate
- fig = plt.figure(figsize=figsize)
- ax = plt.axes()
- clim = [var.min(),var.max()]
- plt.colorbar(mesh.slicer(var, imageType=imageType, normal=normal, index=0, ax=ax, clim=clim))
- tlt = plt.title(normal)
-
- def animateFrame(i):
- mesh.slicer(var, imageType=imageType, normal=normal, index=i, ax=ax, clim=clim)
- tlt.set_text(normal.upper()+('-Slice: %d, %4.4f' % (i,getattr(mesh,'vectorCC'+normal)[i])))
-
- return animate(fig, animateFrame, frames=mesh.vnC['xyz'.index(normal)])
-
- def video(mesh, var, function, figsize=(10, 8), colorbar=True, skip=1):
- """
- Call a function for a list of models to create a video.
-
- ::
-
- def function(var, ax, clim, tlt, i):
- tlt.set_text('%d'%i)
- return mesh.plotImage(var, imageType='CC', ax=ax, clim=clim)
-
- mesh.video([model1, model2, ..., modeln],function)
- """
- # First set up the figure, the axis, and the plot element we want to animate
- fig = plt.figure(figsize=figsize)
- ax = plt.axes()
- VAR = np.concatenate(var)
- clim = [VAR.min(),VAR.max()]
- tlt = plt.title('')
- if colorbar:
- plt.colorbar(function(var[0],ax,clim,tlt,0))
-
- frames = np.arange(0,len(var),skip)
- def animateFrame(j):
- i = frames[j]
- function(var[i],ax,clim,tlt,i)
-
- return animate(fig, animateFrame, frames=len(frames))
class CylView(object):
diff --git a/SimPEG/Tests.py b/SimPEG/Tests.py
index 6d414441..804ca6b8 100644
--- a/SimPEG/Tests.py
+++ b/SimPEG/Tests.py
@@ -1,5 +1,4 @@
import numpy as np
-import matplotlib.pyplot as plt
from numpy.linalg import norm
from SimPEG.Utils import mkvc, sdiag, diagEst
from SimPEG import Utils
@@ -311,6 +310,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
if plotIt:
+ import matplotlib.pyplot as plt
ax = ax or plt.subplot(111)
ax.loglog(h, E0, 'b')
ax.loglog(h, E1, 'g--')
diff --git a/SimPEG/Utils/SolverUtils.py b/SimPEG/Utils/SolverUtils.py
index 279d2b06..26ff3e2a 100644
--- a/SimPEG/Utils/SolverUtils.py
+++ b/SimPEG/Utils/SolverUtils.py
@@ -37,12 +37,20 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
if len(b.shape) == 1 or b.shape[1] == 1:
b = b.flatten()
# Just one RHS
+
+ if b.dtype is np.dtype('O'):
+ b = b.astype(type(b[0]))
+
if factorize:
X = self.solver.solve(b, **self.kwargs)
else:
X = fun(self.A, b, **self.kwargs)
else: # Multiple RHSs
+ if b.dtype is np.dtype('O'):
+ b = b.astype(type(b[0,0]))
+
X = np.empty_like(b)
+
for i in range(b.shape[1]):
if factorize:
X[:,i] = self.solver.solve(b[:,i])
diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py
index 5280ae79..250146c1 100644
--- a/SimPEG/Utils/__init__.py
+++ b/SimPEG/Utils/__init__.py
@@ -3,7 +3,6 @@ from codeutils import *
from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile
from curvutils import volTetra, faceInfo, indexCube
from interputils import interpmat
-from ipythonutils import easyAnimate as animate
from CounterUtils import *
import ModelBuilder
import SolverUtils
diff --git a/SimPEG/Utils/ipythonutils.py b/SimPEG/Utils/ipythonutils.py
deleted file mode 100644
index 7b66e131..00000000
--- a/SimPEG/Utils/ipythonutils.py
+++ /dev/null
@@ -1,28 +0,0 @@
-from tempfile import NamedTemporaryFile
-import matplotlib.pyplot as plt
-from matplotlib import animation
-
-# http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/
-# http://www.renevolution.com/how-to-install-ffmpeg-on-mac-os-x/
-
-VIDEO_TAG = """"""
-
-def anim_to_html(anim):
- if not hasattr(anim, '_encoded_video'):
- with NamedTemporaryFile(suffix='.mp4') as f:
- anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p'])
- video = open(f.name, "rb").read()
- anim._encoded_video = video.encode("base64")
-
- return VIDEO_TAG.format(anim._encoded_video)
-
-def display_animation(anim):
- plt.close(anim._fig)
- return anim_to_html(anim)
-
-animation.Animation._repr_html_ = display_animation
-
-easyAnimate = animation.FuncAnimation
diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py
index ee58bf86..b38bb4a1 100644
--- a/SimPEG/Utils/matutils.py
+++ b/SimPEG/Utils/matutils.py
@@ -399,8 +399,10 @@ def diagEst(matFun, n, k=None, approach='Probing'):
class Zero(object):
def __add__(self, v):return v
def __radd__(self, v):return v
+ def __iadd__(self, v):return v
def __sub__(self, v):return -v
def __rsub__(self, v):return v
+ def __isub__(self, v):return v
def __mul__(self, v):return self
def __rmul__(self, v):return self
def __div__(self, v): return self
diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py
index 369134e7..cc51fd1f 100644
--- a/SimPEG/__init__.py
+++ b/SimPEG/__init__.py
@@ -15,7 +15,7 @@ import Directives
import Inversion
import Tests
-__version__ = '0.1.3'
+__version__ = '0.1.9'
__author__ = 'Rowan Cockett'
__license__ = 'MIT'
__copyright__ = 'Copyright 2014 Rowan Cockett'
diff --git a/docs/api_Examples.rst b/docs/api_Examples.rst
index 214dd8e1..68589fb3 100644
--- a/docs/api_Examples.rst
+++ b/docs/api_Examples.rst
@@ -3,8 +3,15 @@
Examples
********
-Forward problem
-===============
+.. toctree::
+ :maxdepth: 1
+ :glob:
+
+ examples/*
+
+
+External Notebooks
+==================
* `Example 1: Direct Current `_
* `Example 2: Seismic-Acoustic `_
diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst
index 7bf398b1..a7f7abae 100644
--- a/docs/api_Mesh.rst
+++ b/docs/api_Mesh.rst
@@ -23,23 +23,9 @@ the implementations.
.. plot::
- from SimPEG import Mesh, Utils, np
- import matplotlib.pyplot as plt
- sz = [10,10]
- tM = Mesh.TensorMesh(sz)
- qM = Mesh.TreeMesh(sz)
- qM.refine(lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.3 else 0)
- rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
+ from SimPEG import Examples
+ Examples.Mesh_ThreeMeshes.run()
- fig, axes = plt.subplots(1,3,figsize=(14,5))
- opts = {}
- tM.plotGrid(ax=axes[0], **opts)
- axes[0].set_title('TensorMesh')
- qM.plotGrid(ax=axes[1], **opts)
- axes[1].set_title('TreeMesh')
- rM.plotGrid(ax=axes[2], **opts)
- axes[2].set_title('CurvilinearMesh')
- plt.show()
Variable Locations and Terminology
diff --git a/docs/api_Tests.rst b/docs/api_Tests.rst
index 614d8747..b0b2fbd0 100644
--- a/docs/api_Tests.rst
+++ b/docs/api_Tests.rst
@@ -3,6 +3,6 @@
Testing SimPEG
==============
-.. automodule:: SimPEG.Tests.TestUtils
+.. automodule:: SimPEG.Tests
:members:
:undoc-members:
diff --git a/docs/conf.py b/docs/conf.py
index 3bb33621..fee262de 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers'
# built documents.
#
# The short X.Y version.
-version = '0.1.3'
+version = '0.1.9'
# The full version, including alpha/beta/rc tags.
-release = '0.1.3'
+release = '0.1.9'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
diff --git a/docs/examples/EM_FDEM_1D_Inversion.rst b/docs/examples/EM_FDEM_1D_Inversion.rst
new file mode 100644
index 00000000..acbc8cdc
--- /dev/null
+++ b/docs/examples/EM_FDEM_1D_Inversion.rst
@@ -0,0 +1,26 @@
+.. _examples_EM_FDEM_1D_Inversion:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+
+EM: FDEM: 1D: Inversion
+=======================
+
+Here we will create and run a FDEM 1D inversion.
+
+
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.EM_FDEM_1D_Inversion.run()
+
+.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py
+ :language: python
+ :linenos:
diff --git a/docs/examples/FLOW_Richards_1D_Celia1990.rst b/docs/examples/FLOW_Richards_1D_Celia1990.rst
new file mode 100644
index 00000000..d2e01c13
--- /dev/null
+++ b/docs/examples/FLOW_Richards_1D_Celia1990.rst
@@ -0,0 +1,52 @@
+.. _examples_FLOW_Richards_1D_Celia1990:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+
+FLOW: Richards: 1D: Celia1990
+=============================
+
+There are two different forms of Richards equation that differ
+on how they deal with the non-linearity in the time-stepping term.
+
+The most fundamental form, referred to as the
+'mixed'-form of Richards Equation Celia1990_
+
+.. math::
+
+ \frac{\partial \theta(\psi)}{\partial t} - \nabla \cdot k(\psi) \nabla \psi - \frac{\partial k(\psi)}{\partial z} = 0
+ \quad \psi \in \Omega
+
+where \\(\\theta\\) is water content, and \\(\\psi\\) is pressure head.
+This formulation of Richards equation is called the
+'mixed'-form because the equation is parameterized in \\(\\psi\\)
+but the time-stepping is in terms of \\(\\theta\\).
+
+As noted in Celia1990_ the 'head'-based form of Richards
+equation can be written in the continuous form as:
+
+.. math::
+
+ \frac{\partial \theta}{\partial \psi}\frac{\partial \psi}{\partial t} - \nabla \cdot k(\psi) \nabla \psi - \frac{\partial k(\psi)}{\partial z} = 0 \quad \psi \in \Omega
+
+However, it can be shown that this does not conserve mass in the discrete formulation.
+
+Here we reproduce the results from Celia1990_ demonstrating the head-based formulation and the mixed-formulation.
+
+.. _Celia1990: http://www.webpages.uidaho.edu/ch/papers/Celia.pdf
+
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.FLOW_Richards_1D_Celia1990.run()
+
+.. literalinclude:: ../../SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
+ :language: python
+ :linenos:
diff --git a/docs/examples/Forward_BasicDirectCurrent.rst b/docs/examples/Forward_BasicDirectCurrent.rst
new file mode 100644
index 00000000..20b39eb8
--- /dev/null
+++ b/docs/examples/Forward_BasicDirectCurrent.rst
@@ -0,0 +1,21 @@
+.. _examples_Forward_BasicDirectCurrent:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+Forward BasicDirectCurrent
+==========================
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.Forward_BasicDirectCurrent.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Forward_BasicDirectCurrent.py
+ :language: python
+ :linenos:
diff --git a/docs/examples/Inversion_Linear.rst b/docs/examples/Inversion_Linear.rst
new file mode 100644
index 00000000..d635d8e1
--- /dev/null
+++ b/docs/examples/Inversion_Linear.rst
@@ -0,0 +1,26 @@
+.. _examples_Inversion_Linear:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+
+Inversion: Linear Problem
+=========================
+
+Here we go over the basics of creating a linear problem and inversion.
+
+
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.Inversion_Linear.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Inversion_Linear.py
+ :language: python
+ :linenos:
diff --git a/docs/examples/Mesh_Basic_PlotImage.rst b/docs/examples/Mesh_Basic_PlotImage.rst
new file mode 100644
index 00000000..a730f303
--- /dev/null
+++ b/docs/examples/Mesh_Basic_PlotImage.rst
@@ -0,0 +1,27 @@
+.. _examples_Mesh_Basic_PlotImage:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+
+Mesh: Basic: PlotImage
+======================
+
+You can use M.PlotImage to plot images on all of the Meshes.
+
+
+
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.Mesh_Basic_PlotImage.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_PlotImage.py
+ :language: python
+ :linenos:
diff --git a/docs/examples/Mesh_Basic_Types.rst b/docs/examples/Mesh_Basic_Types.rst
new file mode 100644
index 00000000..9bbce0e8
--- /dev/null
+++ b/docs/examples/Mesh_Basic_Types.rst
@@ -0,0 +1,26 @@
+.. _examples_Mesh_Basic_Types:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+
+Mesh: Basic: Types
+==================
+
+Here we show SimPEG used to create three different types of meshes.
+
+
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.Mesh_Basic_Types.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_Types.py
+ :language: python
+ :linenos:
diff --git a/docs/examples/Mesh_Operators_CahnHilliard.rst b/docs/examples/Mesh_Operators_CahnHilliard.rst
new file mode 100644
index 00000000..9786e911
--- /dev/null
+++ b/docs/examples/Mesh_Operators_CahnHilliard.rst
@@ -0,0 +1,57 @@
+.. _examples_Mesh_Operators_CahnHilliard:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+
+Mesh: Operators: Cahn Hilliard
+==============================
+
+This example is based on the example in the FiPy_ library.
+Please see their documentation for more information about the Cahn-Hilliard equation.
+
+The "Cahn-Hilliard" equation separates a field \\( \\phi \\) into 0 and 1 with smooth transitions.
+
+.. math::
+
+ \frac{\partial \phi}{\partial t} = \nabla \cdot D \nabla \left( \frac{\partial f}{\partial \phi} - \epsilon^2 \nabla^2 \phi \right)
+
+Where \\( f \\) is the energy function \\( f = ( a^2 / 2 )\\phi^2(1 - \\phi)^2 \\)
+which drives \\( \\phi \\) towards either 0 or 1, this competes with the term
+\\(\\epsilon^2 \\nabla^2 \\phi \\) which is a diffusion term that creates smooth changes in \\( \\phi \\).
+The equation can be factored:
+
+.. math::
+
+ \frac{\partial \phi}{\partial t} = \nabla \cdot D \nabla \psi \\
+ \psi = \frac{\partial^2 f}{\partial \phi^2} (\phi - \phi^{\text{old}}) + \frac{\partial f}{\partial \phi} - \epsilon^2 \nabla^2 \phi
+
+Here we will need the derivatives of \\( f \\):
+
+.. math::
+
+ \frac{\partial f}{\partial \phi} = (a^2/2)2\phi(1-\phi)(1-2\phi)
+ \frac{\partial^2 f}{\partial \phi^2} = (a^2/2)2[1-6\phi(1-\phi)]
+
+The implementation below uses backwards Euler in time with an exponentially increasing time step.
+The initial \\( \\phi \\) is a normally distributed field with a standard deviation of 0.1 and mean of 0.5.
+The grid is 60x60 and takes a few seconds to solve ~130 times. The results are seen below, and you can see the
+field separating as the time increases.
+
+.. _FiPy: http://www.ctcms.nist.gov/fipy/examples/cahnHilliard/generated/examples.cahnHilliard.mesh2DCoupled.html
+
+
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.Mesh_Operators_CahnHilliard.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_Operators_CahnHilliard.py
+ :language: python
+ :linenos:
diff --git a/docs/examples/Mesh_QuadTree_Creation.rst b/docs/examples/Mesh_QuadTree_Creation.rst
new file mode 100644
index 00000000..5db5a982
--- /dev/null
+++ b/docs/examples/Mesh_QuadTree_Creation.rst
@@ -0,0 +1,31 @@
+.. _examples_Mesh_QuadTree_Creation:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+
+Mesh: QuadTree: Creation
+========================
+
+You can give the refine method a function, which is evaluated on every cell
+of the TreeMesh.
+
+Occasionally it is useful to initially refine to a constant level
+(e.g. 3 in this 32x32 mesh). This means the function is first evaluated
+on an 8x8 mesh (2^3).
+
+
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.Mesh_QuadTree_Creation.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_Creation.py
+ :language: python
+ :linenos:
diff --git a/docs/examples/Mesh_QuadTree_FaceDiv.rst b/docs/examples/Mesh_QuadTree_FaceDiv.rst
new file mode 100644
index 00000000..6bfdd47f
--- /dev/null
+++ b/docs/examples/Mesh_QuadTree_FaceDiv.rst
@@ -0,0 +1,26 @@
+.. _examples_Mesh_QuadTree_FaceDiv:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+
+Mesh: QuadTree: FaceDiv
+=======================
+
+
+
+
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.Mesh_QuadTree_FaceDiv.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_FaceDiv.py
+ :language: python
+ :linenos:
diff --git a/docs/examples/Mesh_QuadTree_HangingNodes.rst b/docs/examples/Mesh_QuadTree_HangingNodes.rst
new file mode 100644
index 00000000..93875478
--- /dev/null
+++ b/docs/examples/Mesh_QuadTree_HangingNodes.rst
@@ -0,0 +1,31 @@
+.. _examples_Mesh_QuadTree_HangingNodes:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+
+Mesh: QuadTree: Hanging Nodes
+=============================
+
+You can give the refine method a function, which is evaluated on every cell
+of the TreeMesh.
+
+Occasionally it is useful to initially refine to a constant level
+(e.g. 3 in this 32x32 mesh). This means the function is first evaluated
+on an 8x8 mesh (2^3).
+
+
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.Mesh_QuadTree_HangingNodes.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_QuadTree_HangingNodes.py
+ :language: python
+ :linenos:
diff --git a/docs/examples/Mesh_Tensor_Creation.rst b/docs/examples/Mesh_Tensor_Creation.rst
new file mode 100644
index 00000000..e6cc5b67
--- /dev/null
+++ b/docs/examples/Mesh_Tensor_Creation.rst
@@ -0,0 +1,43 @@
+.. _examples_Mesh_Tensor_Creation:
+
+.. --------------------------------- ..
+.. ..
+.. THIS FILE IS AUTO GENEREATED ..
+.. ..
+.. SimPEG/Examples/__init__.py ..
+.. ..
+.. --------------------------------- ..
+
+
+
+Mesh: Tensor: Creation
+======================
+
+For tensor meshes, there are some functions that can come
+in handy. For example, creating mesh tensors can be a bit time
+consuming, these can be created speedily by just giving numbers
+and sizes of padding. See the example below, that follows this
+notation::
+
+ h1 = (
+ (cellSize, numPad, [, increaseFactor]),
+ (cellSize, numCore),
+ (cellSize, numPad, [, increaseFactor])
+ )
+
+.. note::
+
+ You can center your mesh by passing a 'C' for the x0[i] position.
+ A 'N' will make the entire mesh negative, and a '0' (or a 0) will
+ make the mesh start at zero.
+
+
+
+.. plot::
+
+ from SimPEG import Examples
+ Examples.Mesh_Tensor_Creation.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_Tensor_Creation.py
+ :language: python
+ :linenos:
diff --git a/setup.py b/setup.py
index 531dc539..bcb5b8e3 100644
--- a/setup.py
+++ b/setup.py
@@ -64,6 +64,7 @@ cython_files = [
"SimPEG/Mesh/TreeUtils"
]
extensions = [Extension(f, [f+ext]) for f in cython_files]
+scripts = [f+'.pyx' for f in cython_files]
if USE_CYTHON and "cleanall" not in args:
from Cython.Build import cythonize
@@ -76,7 +77,7 @@ with open("README.rst") as f:
setup(
name = "SimPEG",
- version = "0.1.3",
+ version = "0.1.9",
packages = find_packages(),
install_requires = ['numpy>=1.7',
'scipy>=0.13',
@@ -95,5 +96,6 @@ setup(
use_2to3 = False,
include_dirs=[np.get_include()],
ext_modules = extensions,
+ scripts=scripts,
**cythonKwargs
)
diff --git a/tests/em/examples/__init__.py b/tests/em/examples/__init__.py
deleted file mode 100644
index 38d84328..00000000
--- a/tests/em/examples/__init__.py
+++ /dev/null
@@ -1,11 +0,0 @@
-if __name__ == '__main__':
- import os
- import glob
- import unittest
- test_file_strings = glob.glob('test_*.py')
- module_strings = [str[0:len(str)-3] for str in test_file_strings]
- suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
- in module_strings]
- testSuite = unittest.TestSuite(suites)
-
- unittest.TextTestRunner(verbosity=2).run(testSuite)
diff --git a/tests/em/examples/test_Examples.py b/tests/em/examples/test_Examples.py
deleted file mode 100644
index 5a601d3b..00000000
--- a/tests/em/examples/test_Examples.py
+++ /dev/null
@@ -1,10 +0,0 @@
-import unittest, os
-from SimPEG.EM import Examples
-
-class EM_ExamplesRunning(unittest.TestCase):
-
- def test_CylInversion(self):
- Examples.CylInversion.run(plotIt=False)
-
-if __name__ == '__main__':
- unittest.main()
diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py
index 52108c4e..d3bcb218 100644
--- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py
+++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py
@@ -38,7 +38,7 @@ def derivTest(fdemType, comp):
survey = prb.survey
def fun(x):
return survey.dpred(x), lambda x: prb.Jvec(x0, x)
- return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
+ return Tests.checkDerivative(fun, x0, num=2, plotIt=False, eps=FLR)
class FDEM_DerivTests(unittest.TestCase):
diff --git a/tests/examples/test_examples.py b/tests/examples/test_examples.py
index 1fcf05f5..2e4803b1 100644
--- a/tests/examples/test_examples.py
+++ b/tests/examples/test_examples.py
@@ -1,17 +1,20 @@
import unittest
import sys
-from SimPEG.Examples import Linear, DCfwd
+from SimPEG import Examples
import numpy as np
-class TestLinear(unittest.TestCase):
- def test_running(self):
- Linear.run(100, plotIt=False)
+def get(test):
+ def test_func(self):
+ print '\nTesting %s.run(plotIt=False)\n'%test
+ getattr(Examples, test).run(plotIt=False)
self.assertTrue(True)
+ return test_func
+attrs = dict()
+for test in Examples.__examples__:
+ attrs['test_'+test] = get(test)
+
+TestExamples = type('TestExamples', (unittest.TestCase,), attrs)
-class TestDCfwd(unittest.TestCase):
- def test_running(self):
- DCfwd.run(plotIt=False)
- self.assertTrue(True)
if __name__ == '__main__':
unittest.main()
diff --git a/tests/flow/test_Richards.py b/tests/flow/test_Richards.py
index 6579370d..d63a6210 100644
--- a/tests/flow/test_Richards.py
+++ b/tests/flow/test_Richards.py
@@ -12,6 +12,8 @@ except Exception, e:
TOL = 1E-8
+np.random.seed(0)
+
class TestModels(unittest.TestCase):
def test_BaseHaverkamp_Theta(self):
diff --git a/tests/flow/test_examples.py b/tests/flow/test_examples.py
deleted file mode 100644
index bc519e6d..00000000
--- a/tests/flow/test_examples.py
+++ /dev/null
@@ -1,12 +0,0 @@
-import unittest
-import sys
-from SimPEG.FLOW.Examples import Celia1990
-import numpy as np
-
-class TestCelia1990(unittest.TestCase):
- def test_running(self):
- Celia1990.run(plotIt=False)
- self.assertTrue(True)
-
-if __name__ == '__main__':
- unittest.main()
diff --git a/tests/utils/test_Zero.py b/tests/utils/test_Zero.py
index be7c9bbe..594de6a6 100644
--- a/tests/utils/test_Zero.py
+++ b/tests/utils/test_Zero.py
@@ -20,6 +20,13 @@ class Tests(unittest.TestCase):
assert 3*z == 0
assert z*3 == 0
assert z/3 == 0
+
+ a = 1
+ a += z
+ assert a == 1
+ a = 1
+ a += z
+ assert a == 1
self.assertRaises(ZeroDivisionError, lambda:3/z)
def test_mat_zero(self):