Added writeUBC_DCobs to StaticUtils from DCIPUtils

Conflicts:
	SimPEG/EM/Base.py
	SimPEG/EM/FDEM/SurveyFDEM.py

SimPEG/Directives.py

@@ -146,10 +146,15 @@ class BetaSchedule(InversionDirective):
 
 class TargetMisfit(InversionDirective):
 
+    chifact = 1.
+    phi_d_star = None
+
     @property
     def target(self):
         if getattr(self, '_target', None) is None:
-            self._target = self.survey.nD*0.5
+            if self.phi_d_star is None:
+                self.phi_d_star = 0.5 * self.survey.nD
+            self._target = self.chifact * self.phi_d_star # the factor of 0.5 is because we do phid = 0.5*|| dpred - dobs||^2
         return self._target
     @target.setter
     def target(self, val):
@@ -305,15 +310,28 @@ class Update_IRLS(InversionDirective):
         self.reg._W = None
 
 class Update_lin_PreCond(InversionDirective):
+    """
+    Create a Jacobi preconditioner for the linear problem
+    """
+    onlyOnStart=False
 
+    def initialize(self):
+
+        if getattr(self.opt, 'approxHinv', None) is None:
+            # Update the pre-conditioner
+            diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
+            self.opt.approxHinv = PC
 
     def endIter(self):
         # Cool the threshold parameter
+        if self.onlyOnStart==True:
+            return
 
         if getattr(self.opt, 'approxHinv', None) is not None:
             # Update the pre-conditioner
             diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
-            PC     = Utils.sdiag(diagA**-1.)
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
             self.opt.approxHinv = PC
 
 

SimPEG/EM/Analytics/DC.py

@@ -0,0 +1,118 @@
+import numpy as np
+from scipy.constants import mu_0, pi
+from scipy import special
+
+def DCAnalyticHalf(txloc, rxlocs, sigma, earth_type="wholespace"):
+    """
+        Analytic solution for electric potential from a postive pole
+
+        :param array txloc: a xyz location of A (+) electrode (np.r_[xa, ya, za])
+        :param list rxlocs: xyz locations of M (+) and N (-) electrodes [M, N]
+
+            e.g.
+            rxlocs = [M, N]
+            M: xyz locations of M (+) electrode (np.c_[xmlocs, ymlocs, zmlocs])
+            N: xyz locations of N (-) electrode (np.c_[xnlocs, ynlocs, znlocs])
+
+        :param float or complex sigma: values of conductivity
+        :param string earth_type: values of conductivity ("wholsespace" or "halfspace")
+
+    """
+    M = rxlocs[0]
+    N = rxlocs[1]
+
+    rM = np.sqrt( (M[:,0]-txloc[0])**2 + (M[:,1]-txloc[1])**2 + (M[:,2]-txloc[1])**2 )
+    rN = np.sqrt( (N[:,0]-txloc[0])**2 + (N[:,1]-txloc[1])**2 + (N[:,2]-txloc[1])**2 )
+
+    phiM = 1./(4*np.pi*rM*sigma)
+    phiN = 1./(4*np.pi*rN*sigma)
+    phi = phiM - phiN
+
+    if earth_type == "halfspace":
+        phi *= 2
+
+    return phi
+
+deg2rad  = lambda deg: deg/180.*np.pi
+rad2deg  = lambda rad: rad*180./np.pi
+
+def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
+                 field_type = "secondary", order=12, halfspace=False):
+# def DCSpherePointCurrent(txloc, rxloc, xc, radius, rho, rho1, \
+#                  field_type = "secondary", order=12):
+    """
+
+        Parameters:
+
+            :param array txloc: A (+) current electrode location (x,y,z)
+            :param array xc: x center of depressed sphere
+            :param array rxloc: M(+) electrode locations / (Nx3 array, # of electrodes)
+
+            :param float radius: radius (float): radius of the sphere (m)
+            :param float rho: resistivity of the background (ohm-m)
+            :param float rho1: resistivity of the sphere
+            :param string field_type: : "secondary", "total", "primary"
+                  (default="secondary")
+                  "secondary": secondary potential only due to sphere
+                  "primary": primary potential from the point source
+                  "total": "secondary"+"primary"
+            :param float order: maximum order of Legendre polynomial (default=12)
+
+        Written by Seogi Kang (skang@eos.ubc.ca)
+        Ph.D. Candidate of University of British Columbia, Canada
+
+    """
+
+    Pleg = []
+    # Compute Legendre Polynomial
+    for i in range(order):
+        Pleg.append(special.legendre(i, monic=0))
+
+
+    rho = 1./sigma
+    rho1 = 1./sigma1
+
+    # Center of the sphere should be aligned in txloc in y-direction
+    yc = txloc[1]
+    xyz = np.c_[rxloc[:,0]-xc, rxloc[:,1]-yc, rxloc[:,2]]
+    r = np.sqrt( (xyz**2).sum(axis=1) )
+
+    x0 = abs(txloc[0]-xc)
+
+    costheta = xyz[:,0]/r * (txloc[0]-xc)/x0
+    phi = np.zeros_like(r)
+    R = (r**2+x0**2.-2.*r*x0*costheta)**0.5
+    # primary potential in a whole space
+    prim = rho*1./(4*np.pi*R)
+
+    if field_type =="primary":
+        return prim
+
+    sphind = r < radius
+    out = np.zeros_like(r)
+    for n in range(order):
+        An, Bn = AnBnfun(n, radius, x0, rho, rho1)
+        dumout = An*r[~sphind]**(-n-1.)*Pleg[n](costheta[~sphind])
+        out[~sphind] += dumout
+        dumin = Bn*r[sphind]**(n)*Pleg[n](costheta[sphind])
+        out[sphind] += dumin
+
+    out[~sphind] += prim[~sphind]
+
+    if halfspace:
+        scale = 2
+    else:
+        scale = 1
+
+    if field_type == "secondary":
+        return scale*(out-prim)
+    elif field_type == "total":
+        return scale*out
+
+def AnBnfun(n, radius, x0, rho, rho1, I=1.):
+    const = I*rho/(4*np.pi)
+    bunmo = n*rho + (n+1)*rho1
+    An = const * radius**(2*n+1) / x0 ** (n+1.) * n * \
+        (rho1-rho) / bunmo
+    Bn = const * 1. / x0 ** (n+1.) * (2*n+1) * (rho1) / bunmo
+    return An, Bn

SimPEG/EM/Analytics/__init__.py

@@ -1,3 +1,4 @@
 from TDEM import hzAnalyticDipoleT
 from FDEM import hzAnalyticDipoleF
 from FDEMcasing import *
+from DC import DCAnalyticHalf, DCAnalyticSphere

SimPEG/EM/Base.py

@@ -1,6 +1,7 @@
 from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
 
+
 class EMPropMap(Maps.PropMap):
     """
         Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m)
@@ -61,6 +62,15 @@ class BaseEMProblem(Problem.BaseProblem):
             self._Me = self.mesh.getEdgeInnerProduct()
         return self._Me
 
+    @property
+    def MeI(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MeI', None) is None:
+            self._MeI = self.mesh.getEdgeInnerProduct(invMat=True)
+        return self._MeI
+
     @property
     def Mf(self):
         """
@@ -70,6 +80,20 @@ class BaseEMProblem(Problem.BaseProblem):
             self._Mf = self.mesh.getFaceInnerProduct()
         return self._Mf
 
+    @property
+    def MfI(self):
+        """
+            Face inner product matrix
+        """
+        if getattr(self, '_MfI', None) is None:
+            self._MfI = self.mesh.getFaceInnerProduct(invMat=True)
+        return self._MfI
+
+    @property
+    def Vol(self):
+        if getattr(self, '_Vol', None) is None:
+            self._Vol = Utils.sdiag(self.mesh.vol)
+        return self._Vol
 
     # ----- Magnetic Permeability ----- #
     @property
@@ -127,7 +151,6 @@ class BaseEMProblem(Problem.BaseProblem):
         """
         return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
 
-
     @property
     def MeSigmaI(self):
         """
@@ -150,7 +173,6 @@ class BaseEMProblem(Problem.BaseProblem):
         return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm))
         # return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u)
 
-
     @property
     def MfRho(self):
         """
@@ -183,7 +205,12 @@ class BaseEMProblem(Problem.BaseProblem):
         """
             Derivative of :code:`MfRhoI` with respect to the model.
         """
-        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
+
+        dMfRhoI_dI = -self.MfRhoI**2
+        dMf_drho = self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u)
+        return dMfRhoI_dI * ( dMf_drho * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) )
+
+        # return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
 
 class BaseEMSurvey(Survey.BaseSurvey):
 
@@ -192,7 +219,7 @@ class BaseEMSurvey(Survey.BaseSurvey):
         self.srcList = srcList
         Survey.BaseSurvey.__init__(self, **kwargs)
 
-    def eval(self, u):
+    def eval(self, f):
         """
         Project fields to receiver locations
         :param Fields u: fields object
@@ -202,8 +229,8 @@ class BaseEMSurvey(Survey.BaseSurvey):
         data = Survey.Data(self)
         for src in self.srcList:
             for rx in src.rxList:
-                data[src, rx] = rx.eval(src, self.mesh, u)
+                data[src, rx] = rx.eval(src, self.mesh, f)
         return data
 
-    def evalDeriv(self, u):
+    def evalDeriv(self, f):
         raise Exception('Use Receivers to project fields deriv.')

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -160,9 +160,9 @@ class Fields(SimPEG.Problem.Fields):
             return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
         return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex)
 
-class Fields_e(Fields):
+class Fields3D_e(Fields):
     """
-    Fields object for Problem_e.
+    Fields object for Problem3D_e.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey
@@ -181,7 +181,7 @@ class Fields_e(Fields):
                   }
 
     def __init__(self, mesh, survey, **kwargs):
-        Fields.__init__(self,mesh,survey,**kwargs)
+        Fields.__init__(self, mesh, survey, **kwargs)
 
     def startup(self):
         self.prob = self.survey.prob
@@ -426,9 +426,9 @@ class Fields_e(Fields):
 
 
 
-class Fields_b(Fields):
+class Fields3D_b(Fields):
     """
-    Fields object for Problem_b.
+    Fields object for Problem3D_b.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey
@@ -693,9 +693,9 @@ class Fields_b(Fields):
         return Zero()
 
 
-class Fields_j(Fields):
+class Fields3D_j(Fields):
     """
-    Fields object for Problem_j.
+    Fields object for Problem3D_j.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey
@@ -988,9 +988,9 @@ class Fields_j(Fields):
         return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
 
 
-class Fields_h(Fields):
+class Fields3D_h(Fields):
     """
-    Fields object for Problem_h.
+    Fields object for Problem3D_h.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey

SimPEG/EM/FDEM/ProblemFDEM.py

@@ -1,7 +1,7 @@
 from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
 from SurveyFDEM import Survey as SurveyFDEM
-from FieldsFDEM import Fields, Fields_e, Fields_b, Fields_h, Fields_j
+from FieldsFDEM import Fields, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
 from SimPEG.EM.Base import BaseEMProblem
 from SimPEG.EM.Utils import omega
 
@@ -17,8 +17,8 @@ class BaseFDEMProblem(BaseEMProblem):
             \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
             {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
 
-        if using the E-B formulation (:code:`Problem_e`
-        or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
+        if using the E-B formulation (:code:`Problem3D_e`
+        or :code:`Problem3D_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
 
         If we write Maxwell's equations in terms of
         \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
@@ -28,7 +28,7 @@ class BaseFDEMProblem(BaseEMProblem):
             \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
             \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
 
-        if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
+        if using the H-J formulation (:code:`Problem3D_j` or :code:`Problem3D_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
 
         The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
     """
@@ -87,7 +87,7 @@ class BaseFDEMProblem(BaseEMProblem):
                 du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
 
                 for rx in src.rxList:
-                    df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                    df_dmFun = getattr(f, '_{0}Deriv'.format(rx.projField), None)
                     df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
                     Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
             Ainv.clean()
@@ -125,7 +125,7 @@ class BaseFDEMProblem(BaseEMProblem):
                 for rx in src.rxList:
                     PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
 
-                    df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                    df_duTFun = getattr(f, '_{0}Deriv'.format(rx.projField), None)
                     df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
 
                     ATinvdf_duT = ATinv * df_duT
@@ -137,10 +137,9 @@ class BaseFDEMProblem(BaseEMProblem):
                     df_dmT = df_dmT + du_dmT
 
                     # TODO: this should be taken care of by the reciever?
-                    real_or_imag = rx.projComp
-                    if real_or_imag is 'real':
+                    if rx.component is 'real':
                         Jtv +=   np.array(df_dmT, dtype=complex).real
-                    elif real_or_imag is 'imag':
+                    elif rx.component is 'imag':
                         Jtv += - np.array(df_dmT, dtype=complex).real
                     else:
                         raise Exception('Must be real or imag')
@@ -167,6 +166,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         for i, src in enumerate(Srcs):
             smi, sei = src.eval(self)
+            #Why are you adding?
             s_m[:,i] = s_m[:,i] + smi
             s_e[:,i] = s_e[:,i] + sei
 
@@ -177,7 +177,7 @@ class BaseFDEMProblem(BaseEMProblem):
 ################################ E-B Formulation #########################################
 ##########################################################################################
 
-class Problem_e(BaseFDEMProblem):
+class Problem3D_e(BaseFDEMProblem):
     """
     By eliminating the magnetic flux density using
 
@@ -199,7 +199,7 @@ class Problem_e(BaseFDEMProblem):
 
     _solutionType = 'eSolution'
     _formulation  = 'EB'
-    fieldsPair    = Fields_e
+    fieldsPair    = Fields3D_e
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -288,7 +288,7 @@ class Problem_e(BaseFDEMProblem):
             return C.T * (MfMui * s_mDeriv(v)) -1j * omega(freq) * s_eDeriv(v)
 
 
-class Problem_b(BaseFDEMProblem):
+class Problem3D_b(BaseFDEMProblem):
     """
     We eliminate :math:`\mathbf{e}` using
 
@@ -310,7 +310,7 @@ class Problem_b(BaseFDEMProblem):
 
     _solutionType = 'bSolution'
     _formulation  = 'EB'
-    fieldsPair    = Fields_b
+    fieldsPair    = Fields3D_b
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -436,7 +436,7 @@ class Problem_b(BaseFDEMProblem):
 ##########################################################################################
 
 
-class Problem_j(BaseFDEMProblem):
+class Problem3D_j(BaseFDEMProblem):
     """
     We eliminate \\\(\\\mathbf{h}\\\) using
 
@@ -458,7 +458,7 @@ class Problem_j(BaseFDEMProblem):
 
     _solutionType = 'jSolution'
     _formulation  = 'HJ'
-    fieldsPair    = Fields_j
+    fieldsPair    = Fields3D_j
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -577,7 +577,7 @@ class Problem_j(BaseFDEMProblem):
 
 
 
-class Problem_h(BaseFDEMProblem):
+class Problem3D_h(BaseFDEMProblem):
     """
     We eliminate \\\(\\\mathbf{j}\\\) using
 
@@ -596,7 +596,7 @@ class Problem_h(BaseFDEMProblem):
 
     _solutionType = 'hSolution'
     _formulation  = 'HJ'
-    fieldsPair    = Fields_h
+    fieldsPair    = Fields3D_h
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)

SimPEG/EM/FDEM/RxFDEM.py

@@ -0,0 +1,126 @@
+import SimPEG
+from SimPEG import sp
+
+class BaseRx(SimPEG.Survey.BaseRx):
+    """
+    Frequency domain receiver base class
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        assert(orientation in ['x','y','z']), "Orientation %s not known. Orientation must be in 'x', 'y', 'z'. Arbitrary orientations have not yet been implemented."%orientation
+        assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not %s"%component
+
+        self.projComp = orientation
+        self.component = component
+
+        SimPEG.Survey.BaseRx.__init__(self, locs, rxType=None) #TODO: remove rxType from baseRx
+
+    def projGLoc(self, u):
+        """Grid Location projection (e.g. Ex Fy ...)"""
+        return u._GLoc(self.projField) + self.projComp
+
+    def eval(self, src, mesh, f):
+        """
+        Project fields to recievers to get data.
+
+        :param Source src: FDEM source
+        :param Mesh mesh: mesh used
+        :param Fields f: fields object
+        :rtype: numpy.ndarray
+        :return: fields projected to recievers
+        """
+
+        P = self.getP(mesh, self.projGLoc(f))
+        f_part_complex = f[src, self.projField]
+        f_part = getattr(f_part_complex, self.component) # get the real or imag component
+
+        return P*f_part
+
+    def evalDeriv(self, src, mesh, f, v, adjoint=False):
+        """
+        Derivative of projected fields with respect to the inversion model times a vector.
+
+        :param Source src: FDEM source
+        :param Mesh mesh: mesh used
+        :param Fields f: fields object
+        :param numpy.ndarray v: vector to multiply
+        :rtype: numpy.ndarray
+        :return: fields projected to recievers
+        """
+
+        P = self.getP(mesh, self.projGLoc(f))
+
+        if not adjoint:
+            Pv_complex = P * v
+            Pv = getattr(Pv_complex, self.component)
+        elif adjoint:
+            Pv_real = P.T * v
+
+            if self.component == 'imag':
+                Pv = 1j*Pv_real
+            elif self.component == 'real':
+                Pv = Pv_real.astype(complex)
+            else:
+                raise NotImplementedError('must be real or imag')
+
+        return Pv
+
+
+class Point_e(BaseRx):
+    """
+    Electric field FDEM receiver
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        self.projField = 'e'
+        super(Point_e, self).__init__(locs, orientation, component)
+
+
+class Point_b(BaseRx):
+    """
+    Magnetic flux FDEM receiver
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        self.projField = 'b'
+        super(Point_b, self).__init__(locs, orientation, component)
+
+
+class Point_h(BaseRx):
+    """
+    Magnetic field FDEM receiver
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        self.projField = 'h'
+        super(Point_h, self).__init__(locs, orientation, component)
+
+
+class Point_j(BaseRx):
+    """
+    Current density FDEM receiver
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        self.projField = 'j'
+        super(Point_j, self).__init__(locs, orientation, component)

SimPEG/EM/FDEM/SrcFDEM.py

@@ -9,12 +9,17 @@ class BaseSrc(Survey.BaseSrc):
     """
 
     freq = None
-    # rxPair = RxFDEM
-    integrate = True
+    integrate = False
+    _ePrimary = None
+    _bPrimary = None
+    _hPrimary = None
+    _jPrimary = None
+
+    def __init__(self, rxList, **kwargs):
+        Survey.BaseSrc.__init__(self, rxList, **kwargs)
 
     def eval(self, prob):
         """
-        Evaluate the source terms.
         - :math:`s_m` : magnetic source term
         - :math:`s_e` : electric source term
 
@@ -51,7 +56,9 @@ class BaseSrc(Survey.BaseSrc):
         :rtype: numpy.ndarray
         :return: primary magnetic flux density
         """
-        return Zero()
+        if self._bPrimary is None:
+            return Zero()
+        return self._bPrimary
 
     def hPrimary(self, prob):
         """
@@ -61,7 +68,9 @@ class BaseSrc(Survey.BaseSrc):
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
-        return Zero()
+        if self._hPrimary is None:
+            return Zero()
+        return self._hPrimary
 
     def ePrimary(self, prob):
         """
@@ -71,7 +80,9 @@ class BaseSrc(Survey.BaseSrc):
         :rtype: numpy.ndarray
         :return: primary electric field
         """
-        return Zero()
+        if self._ePrimary is None:
+            return Zero()
+        return self._ePrimary
 
     def jPrimary(self, prob):
         """
@@ -81,7 +92,9 @@ class BaseSrc(Survey.BaseSrc):
         :rtype: numpy.ndarray
         :return: primary current density
         """
-        return Zero()
+        if self._jPrimary is None:
+            return Zero()
+        return self._jPrimary
 
     def s_m(self, prob):
         """
@@ -136,15 +149,14 @@ class RawVec_e(BaseSrc):
     :param list rxList: receiver list
     :param float freq: frequency
     :param numpy.array s_e: electric source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
 
-    def __init__(self, rxList, freq, s_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
+    def __init__(self, rxList, freq, s_e, **kwargs):
         self._s_e = np.array(s_e, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
 
-        BaseSrc.__init__(self, rxList)
+        BaseSrc.__init__(self, rxList, **kwargs)
 
     def s_e(self, prob):
         """
@@ -166,15 +178,14 @@ class RawVec_m(BaseSrc):
     :param float freq: frequency
     :param rxList: receiver list
     :param numpy.array s_m: magnetic source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
 
-    def __init__(self, rxList, freq, s_m, integrate=True):  #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
+    def __init__(self, rxList, freq, s_m, **kwargs):  #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
         self._s_m = np.array(s_m, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
 
-        BaseSrc.__init__(self, rxList)
+        BaseSrc.__init__(self, rxList, **kwargs)
 
     def s_m(self, prob):
         """
@@ -197,14 +208,13 @@ class RawVec(BaseSrc):
     :param float freq: frequency
     :param numpy.array s_m: magnetic source term
     :param numpy.array s_e: electric source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
-    def __init__(self, rxList, freq, s_m, s_e, integrate=True):
+    def __init__(self, rxList, freq, s_m, s_e, **kwargs):
         self._s_m = np.array(s_m, dtype=complex)
         self._s_e = np.array(s_e, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
-        BaseSrc.__init__(self, rxList)
+        BaseSrc.__init__(self, rxList, **kwargs)
 
     def s_m(self, prob):
         """
@@ -278,14 +288,13 @@ class MagDipole(BaseSrc):
     :param float mu: background magnetic permeability
     """
 
-    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0):
+    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0, **kwargs):
         self.freq = float(freq)
         self.loc = loc
         self.orientation = orientation
         assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
         self.moment = moment
         self.mu = mu
-        self.integrate = False
         BaseSrc.__init__(self, rxList)
 
     def bPrimary(self, prob):
@@ -543,7 +552,7 @@ class CircularLoop(BaseSrc):
             if not prob.mesh.isSymmetric:
                 # TODO ?
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
+            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
             srcfct = MagneticDipoleVectorPotential

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -4,126 +4,9 @@ from SimPEG.EM.Base import BaseEMSurvey
 from scipy.constants import mu_0
 from SimPEG.Utils import Zero, Identity
 import SrcFDEM as Src
+import RxFDEM as Rx
 from SimPEG import sp
 
-
-####################################################
-# Receivers
-####################################################
-
-class Rx(SimPEG.Survey.BaseRx):
-    """
-    Frequency domain receivers
-
-    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
-    :param string rxType: reciever type from knownRxTypes
-    """
-
-    knownRxTypes = {
-                    'exr':['e', 'x', 'real'],
-                    'eyr':['e', 'y', 'real'],
-                    'ezr':['e', 'z', 'real'],
-                    'exi':['e', 'x', 'imag'],
-                    'eyi':['e', 'y', 'imag'],
-                    'ezi':['e', 'z', 'imag'],
-
-                    'bxr':['b', 'x', 'real'],
-                    'byr':['b', 'y', 'real'],
-                    'bzr':['b', 'z', 'real'],
-                    'bxi':['b', 'x', 'imag'],
-                    'byi':['b', 'y', 'imag'],
-                    'bzi':['b', 'z', 'imag'],
-
-                    'jxr':['j', 'x', 'real'],
-                    'jyr':['j', 'y', 'real'],
-                    'jzr':['j', 'z', 'real'],
-                    'jxi':['j', 'x', 'imag'],
-                    'jyi':['j', 'y', 'imag'],
-                    'jzi':['j', 'z', 'imag'],
-
-                    'hxr':['h', 'x', 'real'],
-                    'hyr':['h', 'y', 'real'],
-                    'hzr':['h', 'z', 'real'],
-                    'hxi':['h', 'x', 'imag'],
-                    'hyi':['h', 'y', 'imag'],
-                    'hzi':['h', 'z', 'imag'],
-                   }
-    radius = None
-
-    def __init__(self, locs, rxType):
-        SimPEG.Survey.BaseRx.__init__(self, locs, rxType)
-
-    @property
-    def projField(self):
-        """Field Type projection (e.g. e b ...)"""
-        return self.knownRxTypes[self.rxType][0]
-
-    @property
-    def projComp(self):
-        """Component projection (real/imag)"""
-        return self.knownRxTypes[self.rxType][2]
-
-    def projGLoc(self, u):
-        """Grid Location projection (e.g. Ex Fy ...)"""
-        return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
-
-    def eval(self, src, mesh, f):
-        """
-        Project fields to recievers to get data.
-
-        :param Source src: FDEM source
-        :param Mesh mesh: mesh used
-        :param Fields f: fields object
-        :rtype: numpy.ndarray
-        :return: fields projected to recievers
-        """
-        # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
-        # projGLoc += self.knownRxTypes[self.rxType][1]
-
-        P = self.getP(mesh, self.projGLoc(f))
-        f_part_complex = f[src, self.projField]
-        # get the real or imag component
-        real_or_imag = self.projComp
-        f_part = getattr(f_part_complex, real_or_imag)
-
-        return P*f_part
-
-    def evalDeriv(self, src, mesh, f, v, adjoint=False):
-        """
-        Derivative of projected fields with respect to the inversion model times a vector.
-
-        :param Source src: FDEM source
-        :param Mesh mesh: mesh used
-        :param Fields f: fields object
-        :param numpy.ndarray v: vector to multiply
-        :rtype: numpy.ndarray
-        :return: fields projected to recievers
-        """
-
-        P = self.getP(mesh, self.projGLoc(f))
-
-        if not adjoint:
-            Pv_complex = P * v
-            real_or_imag = self.projComp
-            Pv = getattr(Pv_complex, real_or_imag)
-        elif adjoint:
-            Pv_real = P.T * v
-
-            real_or_imag = self.projComp
-            if real_or_imag == 'imag':
-                Pv = 1j*Pv_real
-            elif real_or_imag == 'real':
-                Pv = Pv_real.astype(complex)
-            else:
-                raise NotImplementedError('must be real or imag')
-
-        return Pv
-
-
-####################################################
-# Survey
-####################################################
-
 class Survey(BaseEMSurvey):
     """
     Frequency domain electromagnetic survey
@@ -132,7 +15,7 @@ class Survey(BaseEMSurvey):
     """
 
     srcPair = Src.BaseSrc
-    rxPair = Rx
+    rxPair = Rx.BaseRx
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency

SimPEG/EM/FDEM/__init__.py

@@ -1,3 +1,5 @@
-from SurveyFDEM import Rx, Src, Survey
-from FDEM import BaseFDEMProblem, Problem_e, Problem_b, Problem_j, Problem_h
-from FieldsFDEM import *
+from SurveyFDEM import Survey
+import SrcFDEM as Src
+import RxFDEM as Rx
+from ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
+from FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h

SimPEG/EM/Static/DC/BoundaryUtils.py

@@ -0,0 +1,160 @@
+import numpy as np
+
+def getxBCyBC_CC(mesh, alpha, beta, gamma):
+# def getxBCyBC(mesh, alpha, beta, gamma):
+    """
+    This is a subfunction generating mixed-boundary condition:
+
+    .. math::
+
+        \nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q
+
+        \rho \vec{j} = -\nabla \phi \phi
+
+        \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega
+
+        xBC = f_1(\alpha, \beta, \gamma)
+        yBC = f(\alpha, \beta, \gamma)
+
+    Computes xBC and yBC for cell-centered discretizations
+    """
+    if mesh.dim == 1: #1D
+        if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2):
+            raise Exception("Lenght of list, alpha should be 2")
+        fCCxm,fCCxp = mesh.cellBoundaryInd
+        nBC = fCCxm.sum()+fCCxp.sum()
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+
+        # h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+        h_xm, h_xp = mesh.hx[0], mesh.hx[-1]
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+
+        xBC = np.r_[xBC_xm, xBC_xp]
+        yBC = np.r_[yBC_xm, yBC_xp]
+
+    elif mesh.dim == 2: #2D
+        if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4):
+            raise Exception("Lenght of list, alpha should be 4")
+
+        fxm,fxp,fym,fyp = mesh.faceBoundaryInd
+        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
+        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
+
+        # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
+        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
+
+        h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
+        h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
+        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
+        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
+        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+        xBC_ym = 0.5*a_ym
+        xBC_yp = 0.5*a_yp/b_yp
+        yBC_ym = 0.5*(1.-b_ym)
+        yBC_yp = 0.5*(1.-1./b_yp)
+
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
+
+        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
+        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
+        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
+        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
+
+        xBC = np.r_[xBC_x, xBC_y]
+        yBC = np.r_[yBC_x, yBC_y]
+
+    elif mesh.dim == 3: #3D
+        if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6):
+            raise Exception("Lenght of list, alpha should be 6")
+        # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd
+        fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd
+        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
+        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
+        alpha_zm, beta_zm, gamma_zm = alpha[4], beta[4], gamma[4]
+        alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5]
+
+        # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
+        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
+        # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2]
+
+        h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
+        h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
+        h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp)
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
+        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
+        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
+        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
+
+        a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm)
+        b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm)
+        a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp)
+        b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+        xBC_ym = 0.5*a_ym
+        xBC_yp = 0.5*a_yp/b_yp
+        yBC_ym = 0.5*(1.-b_ym)
+        yBC_yp = 0.5*(1.-1./b_yp)
+        xBC_zm = 0.5*a_zm
+        xBC_zp = 0.5*a_zp/b_zp
+        yBC_zm = 0.5*(1.-b_zm)
+        yBC_zp = 0.5*(1.-1./b_zp)
+
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
+        sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]])
+
+        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
+        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
+        xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz]
+
+        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
+        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
+        yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz]
+
+        xBC = np.r_[xBC_x, xBC_y, xBC_z]
+        yBC = np.r_[yBC_x, yBC_y, yBC_z]
+
+    return xBC, yBC

SimPEG/EM/Static/DC/FieldsDC.py

@@ -0,0 +1,148 @@
+import SimPEG
+from SimPEG.Utils import Identity, Zero
+import numpy as np
+from scipy.constants import epsilon_0
+
+class Fields(SimPEG.Problem.Fields):
+    knownFields = {}
+    dtype = float
+
+    def _phiDeriv(self, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
+            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
+
+        return np.array(self._phiDeriv_u(src, du_dm_v, adjoint) + self._phiDeriv_m(src, v, adjoint), dtype = float)
+
+    def _eDeriv(self, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
+        return np.array(self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint), dtype = float)
+
+    def _jDeriv(self, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
+        return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = float)
+
+
+class Fields_CC(Fields):
+    knownFields = {'phiSolution':'CC'}
+    aliasFields = {
+                    'phi': ['phiSolution','CC','_phi'],
+                    'j' : ['phiSolution','F','_j'],
+                    'e' : ['phiSolution','F','_e'],
+                    'charge' : ['phiSolution','CC','_charge'],
+                  }
+                  # primary - secondary
+                  # CC variables
+
+    def __init__(self, mesh, survey, **kwargs):
+        Fields.__init__(self, mesh, survey, **kwargs)
+        mesh.setCellGradBC("neumann")
+        cellGrad = mesh.cellGrad
+    def startup(self):
+        self.prob = self.survey.prob
+
+    def _GLoc(self, fieldType):
+        if fieldType == 'phi':
+            return 'CC'
+        elif fieldType == 'e' or fieldType == 'j':
+            return 'F'
+        else:
+            raise Exception('Field type must be phi, e, j')
+
+    def _phi(self, phiSolution, srcList):
+        return phiSolution
+
+    def _phiDeriv_u(self, src, v, adjoint = False):
+        return Identity()*v
+
+    def _phiDeriv_m(self, src, v, adjoint = False):
+        return Zero()
+
+    def _j(self, phiSolution, srcList):
+        """
+            .. math::
+                \mathbf{j} = \mathbf{M}^{f \ -1}_{\rho} \mathbf{G} \phi
+        """
+        return self.prob.MfRhoI*self.prob.Grad*phiSolution
+
+    def _e(self, phiSolution, srcList):
+        """
+            In HJ formulation e is not well-defined!!
+            .. math::
+                \vec{e} = -\nabla \phi
+        """
+        return -self.mesh.cellGrad*phiSolution
+
+    def _charge(self, phiSolution, srcList):
+        """
+            .. math::
+                \int \nabla \codt \vec{e} =  \int \frac{\rho_v }{\epsillon_0}
+        """
+        return epsilon_0*self.prob.Vol*(self.mesh.faceDiv*self._e(phiSolution, srcList))
+
+class Fields_N(Fields):
+    knownFields = {'phiSolution':'N'}
+    aliasFields = {
+                    'phi': ['phiSolution','N','_phi'],
+                    'j' : ['phiSolution','E','_j'],
+                    'e' : ['phiSolution','E','_e'],
+                    'charge' : ['phiSolution','N','_charge'],
+                  }
+                  # primary - secondary
+                  # N variables
+
+    def __init__(self, mesh, survey, **kwargs):
+        Fields.__init__(self, mesh, survey, **kwargs)
+
+    def startup(self):
+        self.prob = self.survey.prob
+
+    def _GLoc(self, fieldType):
+        if fieldType == 'phi':
+            return 'N'
+        elif fieldType == 'e' or fieldType == 'j':
+            return 'E'
+        else:
+            raise Exception('Field type must be phi, e, j')
+
+    def _phi(self, phiSolution, srcList):
+        return phiSolution
+
+    def _phiDeriv_u(self, src, v, adjoint = False):
+        return Identity()*v
+
+    def _phiDeriv_m(self, src, v, adjoint = False):
+        return Zero()
+
+    def _j(self, phiSolution, srcList):
+        """
+            In EB formulation j is not well-defined!!
+            .. math::
+                \mathbf{j} = - \mathbf{M}^{e}_{\sigma} \mathbf{G} \phi
+        """
+        return self.prob.MeSigma * self._e(phiSolution, srcList)
+
+    def _e(self, phiSolution, srcList):
+        """
+            In HJ formulation e is not well-defined!!
+            .. math::
+                \vec{e} = -\nabla \phi
+        """
+        return -self.mesh.nodalGrad * phiSolution
+
+    def _charge(self, phiSolution, srcList):
+        """
+            .. math::
+                \int \nabla \codt \vec{e} =  \int \frac{\rho_v }{\epsillon_0}
+        """
+        return - epsilon_0*(self.mesh.nodalGrad.T*self.mesh.getEdgeInnerProduct()*self._e(phiSolution, srcList))

SimPEG/EM/Static/DC/FieldsDC_2D.py

@@ -0,0 +1,146 @@
+import SimPEG
+from SimPEG.Utils import Identity, Zero
+import numpy as np
+
+class Fields_ky(SimPEG.Problem.TimeFields):
+
+    """
+
+    Fancy Field Storage for a 2.5D code.
+
+    u[:,'phi', kyInd] = phi
+    print u[src0,'phi']
+
+    Only one field type is stored for
+    each problem, the rest are computed. The fields obejct acts like an array and is indexed by
+    .. code-block:: python
+        f = problem.fields(m)
+        e = f[srcList,'e']
+        j = f[srcList,'j']
+
+    If accessing all sources for a given field, use the :code:`:`
+    .. code-block:: python
+        f = problem.fields(m)
+        phi = f[:,'phi']
+        e = f[:,'e']
+        b = f[:,'b']
+    The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies)
+    """
+
+    knownFields = {}
+    dtype = float
+
+    def _phiDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
+            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
+
+        return np.array(self._phiDeriv_u(kyInd, src, du_dm_v, adjoint) + self._phiDeriv_m(kyInd, src, v, adjoint), dtype = float)
+
+    def _eDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
+        return np.array(self._eDeriv_u(kyInd, src, du_dm_v, adjoint) + self._eDeriv_m(kyInd, src, v, adjoint), dtype = float)
+
+    def _jDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)
+        return np.array(self._jDeriv_u(kyInd, src, du_dm_v, adjoint) + self._jDeriv_m(kyInd, src, v, adjoint), dtype = float)
+
+
+    # def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+    #     if adjoint is True:
+    #         return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint)
+    #     return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v)
+
+    # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+    #     if adjoint is True:
+    #         return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
+    #     return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
+
+
+class Fields_ky_CC(Fields_ky):
+    knownFields = {'phiSolution':'CC'}
+    aliasFields = {
+                    'phi': ['phiSolution','CC','_phi'],
+                    'j' : ['phiSolution','F','_j'],
+                    'e' : ['phiSolution','F','_e'],
+                  }
+                  # primary - secondary
+                  # CC variables
+
+    def __init__(self, mesh, survey, **kwargs):
+        Fields_ky.__init__(self, mesh, survey, **kwargs)
+
+    def startup(self):
+        self.prob = self.survey.prob
+
+    def _GLoc(self, fieldType):
+        if fieldType == 'phi':
+            return 'CC'
+        elif fieldType == 'e' or fieldType == 'j':
+            return 'F'
+        else:
+            raise Exception('Field type must be phi, e, j')
+
+    def _phi(self, phiSolution, src, kyInd):
+        return phiSolution
+
+    def _phiDeriv_u(self, kyInd, src, v, adjoint = False):
+        return Identity()*v
+
+    def _phiDeriv_m(self, kyInd, src, v, adjoint = False):
+        return Zero()
+
+    def _j(self, phiSolution, srcList):
+        raise NotImplementedError
+
+    def _e(self, phiSolution, srcList):
+        raise NotImplementedError
+
+class Fields_ky_N(Fields_ky):
+    knownFields = {'phiSolution':'N'}
+    aliasFields = {
+                    'phi': ['phiSolution','N','_phi'],
+                    'j' : ['phiSolution','E','_j'],
+                    'e' : ['phiSolution','E','_e'],
+                  }
+                  # primary - secondary
+                  # CC variables
+
+    def __init__(self, mesh, survey, **kwargs):
+        Fields_ky.__init__(self, mesh, survey, **kwargs)
+
+    def startup(self):
+        self.prob = self.survey.prob
+
+    def _GLoc(self, fieldType):
+        if fieldType == 'phi':
+            return 'N'
+        elif fieldType == 'e' or fieldType == 'j':
+            return 'E'
+        else:
+            raise Exception('Field type must be phi, e, j')
+
+    def _phi(self, phiSolution, src, kyInd):
+        return phiSolution
+
+    def _phiDeriv_u(self, kyInd, src, v, adjoint = False):
+        return Identity()*v
+
+    def _phiDeriv_m(self, kyInd, src, v, adjoint = False):
+        return Zero()
+
+    def _j(self, phiSolution, srcList):
+        raise NotImplementedError
+
+    def _e(self, phiSolution, srcList):
+        raise NotImplementedError

SimPEG/EM/Static/DC/ProblemDC.py

@@ -0,0 +1,298 @@
+from SimPEG import Problem, Utils
+from SimPEG.EM.Base import BaseEMProblem
+from SurveyDC import Survey
+from FieldsDC import Fields, Fields_CC, Fields_N
+from SimPEG.Utils import sdiag
+import numpy as np
+from SimPEG.Utils import Zero
+from BoundaryUtils import getxBCyBC_CC
+
+class BaseDCProblem(BaseEMProblem):
+
+    surveyPair = Survey
+    fieldsPair = Fields
+    Ainv = None
+
+    def fields(self, m):
+        self.curModel = m
+
+        if not self.Ainv == None:
+            self.Ainv.clean()
+
+        f = self.fieldsPair(self.mesh, self.survey)
+        A = self.getA()
+        self.Ainv = self.Solver(A, **self.solverOpts)
+        RHS = self.getRHS()
+        u = self.Ainv * RHS
+        Srcs = self.survey.srcList
+        f[Srcs, self._solutionType] = u
+        return f
+
+    def Jvec(self, m, v, f=None):
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        # Jv = self.dataPair(self.survey) #same size as the data
+        A = self.getA()
+
+        Jv = []
+
+        for src in self.survey.srcList:
+            u_src = f[src, self._solutionType] # solution vector
+            dA_dm_v = self.getADeriv(u_src, v)
+            dRHS_dm_v = self.getRHSDeriv(src, v)
+            du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
+
+            for rx in src.rxList:
+                df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
+                # Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
+                Jv.append(rx.evalDeriv(src, self.mesh, f, df_dm_v))
+        # return Utils.mkvc(Jv)
+        return np.hstack(Jv)
+
+    def Jtvec(self, m, v, f=None):
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        # Ensure v is a data object.
+        if not isinstance(v, self.dataPair):
+            v = self.dataPair(self.survey, v)
+
+        Jtv = np.zeros(m.size)
+        AT = self.getA()
+
+        for src in self.survey.srcList:
+            u_src = f[src, self._solutionType]
+            for rx in src.rxList:
+                PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
+                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
+
+                ATinvdf_duT = self.Ainv * df_duT
+
+                dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
+                dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
+                du_dmT = -dA_dmT + dRHS_dmT
+                Jtv += (df_dmT + du_dmT).astype(float)
+
+        return Utils.mkvc(Jtv)
+
+    def getSourceTerm(self):
+        """
+        takes concept of source and turns it into a matrix
+        """
+        """
+        Evaluates the sources, and puts them in matrix form
+
+        :rtype: (numpy.ndarray, numpy.ndarray)
+        :return: q (nC or nN, nSrc)
+        """
+
+        Srcs = self.survey.srcList
+
+        if self._formulation is 'EB':
+            n = self.mesh.nN
+            # return NotImplementedError
+
+        elif self._formulation is 'HJ':
+            n = self.mesh.nC
+
+        q = np.zeros((n, len(Srcs)))
+
+        for i, src in enumerate(Srcs):
+            q[:,i] = src.eval(self)
+        return q
+
+class Problem3D_CC(BaseDCProblem):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'HJ' # CC potentials means J is on faces
+    fieldsPair    = Fields_CC
+
+    def __init__(self, mesh, **kwargs):
+        BaseDCProblem.__init__(self, mesh, **kwargs)
+        self.setBC()
+
+    def getA(self):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = D MfRhoI G
+
+        """
+
+        D = self.Div
+        G = self.Grad
+        MfRhoI = self.MfRhoI
+        A = D * MfRhoI * G
+
+        # I think we should deprecate this for DC problem.
+        # if self._makeASymmetric is True:
+        #     return V.T * A
+        return A
+
+    def getADeriv(self, u, v, adjoint= False):
+
+        D = self.Div
+        G = self.Grad
+        MfRhoIDeriv = self.MfRhoIDeriv
+
+        if adjoint:
+            return(MfRhoIDeriv( G * u ).T) * ( D.T * v)
+
+        return D * (MfRhoIDeriv( G * u ) * v)
+
+    def getRHS(self):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm()
+
+        return RHS
+
+    def getRHSDeriv(self, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
+        # return qDeriv
+        return Zero()
+
+    def setBC(self):
+        if self.mesh.dim==3:
+            fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+            gBFzm = self.mesh.gridFz[fzm,:]
+            gBFzp = self.mesh.gridFz[fzp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+            temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+            alpha_zm, alpha_zp = temp_zm*0., temp_zp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+            beta_zm, beta_zp = temp_zm, temp_zp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+            gamma_zm, gamma_zp = temp_zm*0., temp_zp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
+
+        elif self.mesh.dim==2:
+
+            fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
+
+        x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma)
+        V = self.Vol
+        self.Div = V * self.mesh.faceDiv
+        P_BC, B = self.mesh.getBCProjWF_simple()
+        M = B*self.mesh.aveCC2F
+        self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
+
+
+class Problem3D_N(BaseDCProblem):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'EB' # N potentials means B is on faces
+    fieldsPair    = Fields_N
+
+    def __init__(self, mesh, **kwargs):
+        BaseDCProblem.__init__(self, mesh, **kwargs)
+
+    def getA(self):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = G.T MeSigma G
+
+        """
+
+        MeSigma = self.MeSigma
+        Grad = self.mesh.nodalGrad
+        A = Grad.T * MeSigma * Grad
+
+        # Handling Null space of A
+        A[0,0] = A[0,0] + 1.
+
+        return A
+
+    def getADeriv(self, u, v, adjoint=False):
+        """
+
+        Product of the derivative of our system matrix with respect to the model and a vector
+
+        """
+        MeSigma = self.MeSigma
+        Grad = self.mesh.nodalGrad
+        if not adjoint:
+            return Grad.T*(self.MeSigmaDeriv(Grad*u)*v)
+        elif adjoint:
+            return self.MeSigmaDeriv(Grad*u).T * (Grad*v)
+
+
+    def getRHS(self):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm()
+        return RHS
+
+    def getRHSDeriv(self, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
+        # return qDeriv
+        return Zero()
+
+
+
+

SimPEG/EM/Static/DC/ProblemDC_2D.py

@@ -0,0 +1,349 @@
+from SimPEG import Problem, Utils
+from SimPEG.EM.Base import BaseEMProblem
+from SurveyDC import Survey, Survey_ky
+from FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
+from SimPEG.Utils import sdiag
+import numpy as np
+from SimPEG.Utils import Zero
+from BoundaryUtils import getxBCyBC_CC
+
+class BaseDCProblem_2D(BaseEMProblem):
+
+    surveyPair = Survey_ky
+    fieldsPair = Fields_ky
+    nky = 15
+    kys = np.logspace(-4, 1, nky)
+    Ainv = [None for i in range(nky)]
+    nT = nky # Only for using TimeFields
+
+    def fields(self, m):
+        self.curModel = m
+
+        if not self.Ainv[0] == None:
+            for i in range(self.nky):
+                self.Ainv[i].clean()
+
+        f = self.fieldsPair(self.mesh, self.survey)
+        Srcs = self.survey.srcList
+        for iky in range(self.nky):
+            ky = self.kys[iky]
+            A = self.getA(ky)
+            self.Ainv[iky] = self.Solver(A, **self.solverOpts)
+            RHS = self.getRHS(ky)
+            u = self.Ainv[iky] * RHS
+            f[Srcs, self._solutionType, iky] = u
+        return f
+
+    def Jvec(self, m, v, f=None):
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        Jv = self.dataPair(self.survey) #same size as the data
+        Jv0 = self.dataPair(self.survey)
+
+        # Assume y=0.
+        # This needs some thoughts to implement in general when src is dipole
+        dky = np.diff(self.kys)
+        dky = np.r_[dky[0], dky]
+        y = 0.
+
+        #TODO: this loop is pretty slow .. (Parellize)
+        for iky in range(self.nky):
+            ky = self.kys[iky]
+            A = self.getA(ky)
+            for src in self.survey.srcList:
+                u_src = f[src, self._solutionType, iky] # solution vector
+                dA_dm_v = self.getADeriv(ky, u_src, v)
+                dRHS_dm_v = self.getRHSDeriv(ky, src, v)
+                du_dm_v = self.Ainv[iky] * ( - dA_dm_v + dRHS_dm_v )
+                for rx in src.rxList:
+                    df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                    df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False)
+                    # Trapezoidal intergration
+                    Jv1_temp = 1./np.pi*rx.evalDeriv(ky, src, self.mesh, f, df_dm_v)
+                    if iky==0:
+                        #First assigment
+                        Jv[src, rx] = Jv1_temp*dky[iky]*np.cos(ky*y)
+                    else:
+                        Jv[src, rx] += Jv1_temp*dky[iky]    /2.*np.cos(ky*y)
+                        Jv[src, rx] += Jv0[src, rx]*dky[iky]/2.*np.cos(ky*y)
+                    Jv0[src, rx] = Jv1_temp.copy()
+        return Utils.mkvc(Jv)
+
+    def Jtvec(self, m, v, f=None):
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        # Ensure v is a data object.
+        if not isinstance(v, self.dataPair):
+            v = self.dataPair(self.survey, v)
+
+        Jtv = np.zeros(m.size, dtype=float)
+
+        # Assume y=0.
+        # This needs some thoughts to implement in general when src is dipole
+        dky = np.diff(self.kys)
+        dky = np.r_[dky[0], dky]
+        y = 0.
+
+        for src in self.survey.srcList:
+            for rx in src.rxList:
+                Jtv_temp1 = np.zeros(m.size, dtype=float)
+                Jtv_temp0 = np.zeros(m.size, dtype=float)
+                #TODO: this loop is pretty slow .. (Parellize)
+                for iky in range(self.nky):
+                    u_src = f[src, self._solutionType, iky]
+                    ky = self.kys[iky]
+                    AT = self.getA(ky)
+                    PTv = rx.evalDeriv(ky, src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
+                    df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                    df_duT, df_dmT = df_duTFun(iky, src, None, PTv, adjoint=True)
+
+                    ATinvdf_duT = self.Ainv[iky] * df_duT
+
+                    dA_dmT = self.getADeriv(ky, u_src, ATinvdf_duT, adjoint=True)
+                    dRHS_dmT = self.getRHSDeriv(ky, src, ATinvdf_duT, adjoint=True)
+                    du_dmT = -dA_dmT + dRHS_dmT
+                    Jtv_temp1 = 1./np.pi*(df_dmT + du_dmT).astype(float)
+                    # Trapezoidal intergration
+                    if iky==0:
+                        #First assigment
+                        Jtv += Jtv_temp1*dky[iky]*np.cos(ky*y)
+                    else:
+                        Jtv += Jtv_temp1*dky[iky]/2.*np.cos(ky*y)
+                        Jtv += Jtv_temp0*dky[iky]/2.*np.cos(ky*y)
+                    Jtv_temp0 = Jtv_temp1.copy()
+        return Utils.mkvc(Jtv)
+
+    def getSourceTerm(self, ky):
+        """
+        takes concept of source and turns it into a matrix
+        """
+        """
+        Evaluates the sources, and puts them in matrix form
+
+        :rtype: (numpy.ndarray, numpy.ndarray)
+        :return: q (nC or nN, nSrc)
+        """
+
+        Srcs = self.survey.srcList
+
+        if self._formulation is 'EB':
+            n = self.mesh.nN
+            # return NotImplementedError
+
+        elif self._formulation is 'HJ':
+            n = self.mesh.nC
+
+        q = np.zeros((n, len(Srcs)))
+
+        for i, src in enumerate(Srcs):
+            q[:,i] = src.eval(self)
+        return q
+
+class Problem2D_CC(BaseDCProblem_2D):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'HJ' # CC potentials means J is on faces
+    fieldsPair    = Fields_ky_CC
+
+    def __init__(self, mesh, **kwargs):
+        BaseDCProblem_2D.__init__(self, mesh, **kwargs)
+        self.setBC()
+
+    def getA(self, ky):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = D MfRhoI G
+
+        """
+
+        D = self.Div
+        G = self.Grad
+        vol = self.mesh.vol
+        MfRhoI = self.MfRhoI
+        # Get resistivity rho
+        rho = self.curModel.rho
+        A = D * MfRhoI * G + Utils.sdiag(ky**2*vol/rho)
+        return A
+
+    def getADeriv(self, ky, u, v, adjoint= False):
+
+        D = self.Div
+        G = self.Grad
+        vol = self.mesh.vol
+        MfRhoIDeriv = self.MfRhoIDeriv
+        rho = self.curModel.rho
+        if adjoint:
+            return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v
+        return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v
+
+    def getRHS(self, ky):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm(ky)
+        return RHS
+
+    def getRHSDeriv(self, ky, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, ky, adjoint=adjoint)
+        # return qDeriv
+        return Zero()
+
+    def setBC(self):
+        if self.mesh.dim==3:
+            fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+            gBFzm = self.mesh.gridFz[fzm,:]
+            gBFzp = self.mesh.gridFz[fzp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+            temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+            alpha_zm, alpha_zp = temp_zm*0., temp_zp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+            beta_zm, beta_zp = temp_zm, temp_zp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+            gamma_zm, gamma_zp = temp_zm*0., temp_zp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
+
+        elif self.mesh.dim==2:
+
+            fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
+
+        x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma)
+        V = self.Vol
+        self.Div = V * self.mesh.faceDiv
+        P_BC, B = self.mesh.getBCProjWF_simple()
+        M = B*self.mesh.aveCC2F
+        self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
+
+class Problem2D_N(BaseDCProblem_2D):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'EB' # CC potentials means J is on faces
+    fieldsPair    = Fields_ky_N
+
+    def __init__(self, mesh, **kwargs):
+        BaseDCProblem_2D.__init__(self, mesh, **kwargs)
+        # self.setBC()
+
+    @property
+    def MnSigma(self):
+        """
+            Node inner product matrix for \\(\\sigma\\). Used in the E-B formulation
+        """
+        # TODO: only works isotropic sigma
+        sigma = self.curModel.sigma
+        vol = self.mesh.vol
+        MnSigma = Utils.sdiag(self.mesh.aveN2CC.T*(Utils.sdiag(vol)*sigma))
+
+        return MnSigma
+
+    def MnSigmaDeriv(self, u):
+        """
+            Derivative of MnSigma with respect to the model
+        """
+        sigma = self.curModel.sigma
+        sigmaderiv = self.curModel.sigmaDeriv
+        vol = self.mesh.vol
+        return Utils.sdiag(u)*self.mesh.aveN2CC.T*Utils.sdiag(vol) * self.curModel.sigmaDeriv
+
+    def getA(self, ky):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = D MfRhoI G
+
+        """
+
+        MeSigma = self.MeSigma
+        MnSigma = self.MnSigma
+        Grad = self.mesh.nodalGrad
+        # Get conductivity sigma
+        sigma = self.curModel.sigma
+        A = Grad.T * MeSigma * Grad + ky**2*MnSigma
+
+        # Handling Null space of A
+        A[0,0] = A[0,0] + 1.
+        return A
+
+    def getADeriv(self, ky, u, v, adjoint= False):
+
+        MeSigma = self.MeSigma
+        Grad = self.mesh.nodalGrad
+        sigma = self.curModel.sigma
+        vol = self.mesh.vol
+
+        if adjoint:
+            return self.MeSigmaDeriv(Grad*u).T * (Grad*v) + ky**2*self.MnSigmaDeriv(u).T*v
+        return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) + ky**2*self.MnSigmaDeriv(u)*v
+
+    def getRHS(self, ky):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm(ky)
+        return RHS
+
+    def getRHSDeriv(self, ky, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, ky, adjoint=adjoint)
+        # return qDeriv
+        return Zero()

SimPEG/EM/Static/DC/RxDC.py

@@ -0,0 +1,129 @@
+import SimPEG
+import numpy as np
+from SimPEG.Utils import Zero, closestPoints
+
+class BaseRx(SimPEG.Survey.BaseRx):
+    locs = None
+    rxType = None
+
+    knownRxTypes = {
+                    'phi':['phi',None],
+                    'ex':['e','x'],
+                    'ey':['e','y'],
+                    'ez':['e','z'],
+                    'jx':['j','x'],
+                    'jy':['j','y'],
+                    'jz':['j','z'],
+                    }
+
+    def __init__(self, locs, rxType, **kwargs):
+        SimPEG.Survey.BaseRx.__init__(self, locs, rxType, **kwargs)
+
+
+    @property
+    def projField(self):
+        """Field Type projection (e.g. e b ...)"""
+        return self.knownRxTypes[self.rxType][0]
+
+    def projGLoc(self, f):
+        """Grid Location projection (e.g. Ex Fy ...)"""
+        comp = self.knownRxTypes[self.rxType][1]
+        if comp is not None:
+            return f._GLoc(self.rxType) + comp
+        return f._GLoc(self.rxType)
+
+    def eval(self, src, mesh, f):
+        P = self.getP(mesh, self.projGLoc(f))
+        return P*f[src, self.projField]
+
+    def evalDeriv(self, src, mesh, f, v, adjoint=False):
+        P = self.getP(mesh, self.projGLoc(f))
+        if not adjoint:
+            return P*v
+        elif adjoint:
+            return P.T*v
+
+# DC.Rx.Dipole(locs)
+class Dipole(BaseRx):
+
+    def __init__(self, locsM, locsN, rxType = 'phi', **kwargs):
+        assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size'
+        locs = [locsM, locsN]
+        # We may not need this ...
+        BaseRx.__init__(self, locs, rxType)
+
+    @property
+    def nD(self):
+        """Number of data in the receiver."""
+        return self.locs[0].shape[0]
+
+        # Not sure why ...
+        # return int(self.locs[0].size / 2)
+
+
+    def getP(self, mesh, Gloc):
+        if mesh in self._Ps:
+            return self._Ps[mesh]
+
+        P0 = mesh.getInterpolationMat(self.locs[0], Gloc)
+        P1 = mesh.getInterpolationMat(self.locs[1], Gloc)
+        P = P0 - P1
+
+        if self.storeProjections:
+            self._Ps[mesh] = P
+
+        return P
+
+
+class Dipole_ky(BaseRx):
+
+    def __init__(self, locsM, locsN, rxType = 'phi', **kwargs):
+        assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size'
+        locs = [locsM, locsN]
+        # We may not need this ...
+        BaseRx.__init__(self, locs, rxType)
+
+    @property
+    def nD(self):
+        """Number of data in the receiver."""
+        return self.locs[0].shape[0]
+
+        # Not sure why ...
+        # return int(self.locs[0].size / 2)
+
+    def getP(self, mesh, Gloc):
+        if mesh in self._Ps:
+            return self._Ps[mesh]
+
+        P0 = mesh.getInterpolationMat(self.locs[0], Gloc)
+        P1 = mesh.getInterpolationMat(self.locs[1], Gloc)
+        P = P0 - P1
+        if self.storeProjections:
+            self._Ps[mesh] = P
+        return P
+
+    def eval(self, kys, src, mesh, f):
+        P = self.getP(mesh, self.projGLoc(f))
+        Pf = P*f[src, self.projField,:]
+        return self.IntTrapezoidal(kys, Pf, y=0.)
+
+    def evalDeriv(self, ky, src, mesh, f, v, adjoint=False):
+        P = self.getP(mesh, self.projGLoc(f))
+        if not adjoint:
+            return P*v
+        elif adjoint:
+            return P.T*v
+
+    def IntTrapezoidal(self, kys, Pf, y=0.):
+        phi = np.zeros(Pf.shape[0])
+        nky = kys.size
+        dky = np.diff(kys)
+        dky = np.r_[dky[0], dky]
+        phi0 = 1./np.pi*Pf[:,0]
+        for iky in range(nky):
+            phi1 = 1./np.pi*Pf[:,iky]
+            phi += phi1*dky[iky]/2.*np.cos(kys[iky]*y)
+            phi += phi0*dky[iky]/2.*np.cos(kys[iky]*y)
+            phi0 = phi1.copy()
+        return phi
+

SimPEG/EM/Static/DC/SrcDC.py

@@ -0,0 +1,86 @@
+import SimPEG
+# from SimPEG.EM.Base import BaseEMSurvey
+from SimPEG.Utils import Zero, closestPoints, mkvc
+import numpy as np
+
+class BaseSrc(SimPEG.Survey.BaseSrc):
+
+    current = 1.0
+    loc = None
+
+    def __init__(self, rxList, **kwargs):
+        SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs)
+
+    def eval(self, prob):
+        raise NotImplementedError
+
+    def evalDeriv(self, prob):
+        return Zero()
+
+
+class Dipole(BaseSrc):
+
+    def __init__(self, rxList, locA, locB, **kwargs):
+        assert locA.shape == locB.shape, 'Shape of locA and locB should be the same'
+        self.loc = [locA, locB]
+        BaseSrc.__init__(self, rxList, **kwargs)
+
+    def eval(self, prob):
+        if prob._formulation == 'HJ':
+            inds = closestPoints(prob.mesh, self.loc, gridLoc='CC')
+            q = np.zeros(prob.mesh.nC)
+            q[inds] = self.current * np.r_[1., -1.]
+        elif prob._formulation == 'EB':
+            qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense()
+            qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense()
+            q = self.current * mkvc(qa+qb)
+        return q
+
+class Pole(BaseSrc):
+
+    def __init__(self, rxList, loc, **kwargs):
+        BaseSrc.__init__(self, rxList, loc=loc, **kwargs)
+
+    def eval(self, prob):
+        if prob._formulation == 'HJ':
+            inds = closestPoints(prob.mesh, self.loc)
+            q = np.zeros(prob.mesh.nC)
+            q[inds] = self.current * np.r_[1.]
+        elif prob._formulation == 'EB':
+            q = prob.mesh.getInterpolationMat(self.loc, locType='N').todense()
+            q = self.current * mkvc(q)
+        return q
+
+
+# class Dipole_ky(BaseSrc):
+
+#     def __init__(self, rxList, locA, locB, **kwargs):
+#         assert locA.shape == locB.shape, 'Shape of locA and locB should be the same'
+#         self.loc = [locA[[0,2]], locB[[0,2]]]
+#         BaseSrc.__init__(self, rxList, **kwargs)
+
+#     def eval(self, prob):
+#         if prob._formulation == 'HJ':
+#             inds = closestPoints(prob.mesh, self.loc, gridLoc='CC')
+#             q = np.zeros(prob.mesh.nC)
+#             q[inds] = self.current * np.r_[1., -1.]
+#         elif prob._formulation == 'EB':
+#             qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense()
+#             qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense()
+#             q = self.current * mkvc(qa+qb)
+#         return q
+
+# class Pole_ky(BaseSrc):
+
+#     def __init__(self, rxList, loc, **kwargs):
+#         BaseSrc.__init__(self, rxList, loc=loc, **kwargs)
+
+#     def eval(self, prob):
+#         if prob._formulation == 'HJ':
+#             inds = closestPoints(prob.mesh, self.loc[[0,2]])
+#             q = np.zeros(prob.mesh.nC)
+#             q[inds] = self.current * np.r_[1.]
+#         elif prob._formulation == 'EB':
+#             q = prob.mesh.getInterpolationMat(self.loc[[0,2]], locType='N').todense()
+#             q = self.current * mkvc(q)
+#         return q

SimPEG/EM/Static/DC/SurveyDC.py

@@ -0,0 +1,38 @@
+import SimPEG
+from SimPEG.EM.Base import BaseEMSurvey
+from SimPEG import sp, Survey
+from SimPEG.Utils import Zero, Identity
+from RxDC import BaseRx
+from SrcDC import BaseSrc
+
+class Survey(BaseEMSurvey):
+    rxPair  = BaseRx
+    srcPair = BaseSrc
+
+    def __init__(self, srcList, **kwargs):
+        self.srcList = srcList
+        BaseEMSurvey.__init__(self, srcList, **kwargs)
+
+class Survey_ky(BaseEMSurvey):
+    rxPair  = BaseRx
+    srcPair = BaseSrc
+
+    def __init__(self, srcList, **kwargs):
+        self.srcList = srcList
+        BaseEMSurvey.__init__(self, srcList, **kwargs)
+
+    def eval(self, f):
+        """
+        Project fields to receiver locations
+        :param Fields u: fields object
+        :rtype: numpy.ndarray
+        :return: data
+        """
+        data = SimPEG.Survey.Data(self)
+        kys = self.prob.kys
+        for src in self.srcList:
+            for rx in src.rxList:
+                data[src, rx] = rx.eval(kys, src, self.mesh, f)
+        return data
+
+

SimPEG/EM/Static/DC/Utils.py

@@ -0,0 +1,38 @@
+import numpy as np
+
+def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
+
+    import SimPEG.EM.Static.DC as DC
+
+    elocs = np.arange(0,aSpacing*nElecs,aSpacing)
+    elocs -= (nElecs*aSpacing - aSpacing)/2
+    space = 1
+    WENNER = np.zeros((0,),dtype=int)
+    for ii in range(nElecs):
+        for jj in range(nElecs):
+            test = np.r_[jj,jj+space,jj+space*2,jj+space*3]
+            if np.any(test >= nElecs):
+                break
+            WENNER = np.r_[WENNER, test]
+        space += 1
+    WENNER = WENNER.reshape((-1,4))
+
+
+    if plotIt:
+        for i, s in enumerate('rbkg'):
+            plt.plot(elocs[WENNER[:,i]],s+'.')
+        plt.show()
+
+    # Create sources and receivers
+    i = 0
+    if in2D:
+        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0]
+    else:
+        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0, 0]
+    srcList = []
+    for i in range(WENNER.shape[0]):
+        rx = DC.Rx.Dipole(getLoc(i,1).reshape([1,-1]),getLoc(i,2).reshape([1,-1]))
+        src = DC.Src.Dipole([rx], getLoc(i,0),getLoc(i,3))
+        srcList += [src]
+
+    return srcList

SimPEG/EM/Static/DC/__init__.py

@@ -0,0 +1,8 @@
+from ProblemDC import Problem3D_CC, Problem3D_N
+from ProblemDC_2D import Problem2D_CC, Problem2D_N
+from SurveyDC import Survey, Survey_ky
+import SrcDC as Src #Pole
+import RxDC as Rx
+from FieldsDC import Fields_CC
+from BoundaryUtils import getxBCyBC_CC
+import Utils

SimPEG/EM/Static/IP/ProblemIP.py

@@ -0,0 +1,376 @@
+from SimPEG import Problem, Utils, Maps, Mesh
+from SimPEG.EM.Base import BaseEMProblem
+from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
+from SimPEG.Utils import sdiag
+import numpy as np
+from SimPEG.Utils import Zero
+from SimPEG.EM.Static.DC import getxBCyBC_CC
+from SurveyIP import Survey
+
+class IPPropMap(Maps.PropMap):
+    """
+        Property Map for IP Problems. The electrical chargeability,
+        (\\(\\eta\\)) is the default inversion property
+    """
+    eta = Maps.Property("Electrical Chargeability", defaultInvProp = True)
+
+class BaseIPProblem(BaseEMProblem):
+
+    surveyPair = Survey
+    fieldsPair = Fields
+    PropMap = IPPropMap
+    Ainv = None
+    sigma = None
+    rho = None
+    f = None
+    Ainv = None
+
+    def fields(self, m):
+        self.curModel = m
+        if self.f is None:
+            self.f = self.fieldsPair(self.mesh, self.survey)
+            if self.Ainv == None:
+                A = self.getA()
+                self.Ainv = self.Solver(A, **self.solverOpts)
+            RHS = self.getRHS()
+            u = self.Ainv * RHS
+            Srcs = self.survey.srcList
+            self.f[Srcs, self._solutionType] = u
+        return self.f
+
+    def Jvec(self, m, v, f=None):
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        # Jv = self.dataPair(self.survey) #same size as the data
+        Jv = []
+
+        A = self.getA()
+
+        for src in self.survey.srcList:
+            u_src = f[src, self._solutionType] # solution vector
+            dA_dm_v = self.getADeriv(u_src, v)
+            dRHS_dm_v = self.getRHSDeriv(src, v)
+            du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
+
+            for rx in src.rxList:
+                df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
+                # Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
+                Jv.append(rx.evalDeriv(src, self.mesh, f, df_dm_v))
+        # Conductivity (d u / d log sigma)
+        if self._formulation is 'EB':
+            # return -Utils.mkvc(Jv)
+            return -np.hstack(Jv)
+        # Conductivity (d u / d log rho)
+        if self._formulation is 'HJ':
+            # return Utils.mkvc(Jv)
+            return np.hstack(Jv)
+
+    def Jtvec(self, m, v, f=None):
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        # Ensure v is a data object.
+        if not isinstance(v, self.dataPair):
+            v = self.dataPair(self.survey, v)
+
+        Jtv = np.zeros(m.size)
+        AT = self.getA()
+
+        for src in self.survey.srcList:
+            u_src = f[src, self._solutionType]
+            for rx in src.rxList:
+                PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
+                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
+                ATinvdf_duT = self.Ainv * df_duT
+                dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
+                dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
+                du_dmT = -dA_dmT + dRHS_dmT
+                Jtv += (df_dmT + du_dmT).astype(float)
+        # Conductivity ((d u / d log sigma).T)
+        if self._formulation is 'EB':
+            return -Utils.mkvc(Jtv)
+        # Conductivity ((d u / d log rho).T)
+        if self._formulation is 'HJ':
+            return Utils.mkvc(Jtv)
+
+    def getSourceTerm(self):
+        """
+        takes concept of source and turns it into a matrix
+        """
+        """
+        Evaluates the sources, and puts them in matrix form
+
+        :rtype: (numpy.ndarray, numpy.ndarray)
+        :return: q (nC or nN, nSrc)
+        """
+
+        Srcs = self.survey.srcList
+
+        if self._formulation is 'EB':
+            n = self.mesh.nN
+            # return NotImplementedError
+
+        elif self._formulation is 'HJ':
+            n = self.mesh.nC
+
+        q = np.zeros((n, len(Srcs)))
+
+        for i, src in enumerate(Srcs):
+            q[:,i] = src.eval(self)
+        return q
+
+    @property
+    def deleteTheseOnModelUpdate(self):
+        toDelete = []
+        return toDelete
+
+    # assume log rho or log cond
+    @property
+    def MeSigma(self):
+        """
+            Edge inner product matrix for \\(\\sigma\\). Used in the E-B formulation
+        """
+        if getattr(self, '_MeSigma', None) is None:
+            self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma)
+        return self._MeSigma
+
+    @property
+    def MfRhoI(self):
+        """
+            Inverse of :code:`MfRho`
+        """
+        if getattr(self, '_MfRhoI', None) is None:
+            self._MfRhoI = self.mesh.getFaceInnerProduct(self.rho, invMat=True)
+        return self._MfRhoI
+
+    def MfRhoIDeriv(self,u):
+        """
+            Derivative of :code:`MfRhoI` with respect to the model.
+        """
+
+        dMfRhoI_dI = -self.MfRhoI**2
+        dMf_drho = self.mesh.getFaceInnerProductDeriv(self.rho)(u)
+        drho_dlogrho = Utils.sdiag(self.rho)*self.curModel.etaDeriv
+        return dMfRhoI_dI * ( dMf_drho * ( drho_dlogrho))
+
+    # TODO: This should take a vector
+    def MeSigmaDeriv(self, u):
+        """
+            Derivative of MeSigma with respect to the model
+        """
+        dsigma_dlogsigma = Utils.sdiag(self.sigma)*self.curModel.etaDeriv
+        return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) * dsigma_dlogsigma
+
+class Problem3D_CC(BaseIPProblem):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'HJ' # CC potentials means J is on faces
+    fieldsPair    = Fields_CC
+
+    def __init__(self, mesh, **kwargs):
+        BaseIPProblem.__init__(self, mesh, **kwargs)
+        self.setBC()
+
+    def getA(self):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = D MfRhoI G
+
+        """
+
+        D = self.Div
+        G = self.Grad
+        MfRhoI = self.MfRhoI
+        A = D * MfRhoI * G
+
+        # I think we should deprecate this for DC problem.
+        # if self._makeASymmetric is True:
+        #     return V.T * A
+        return A
+
+    def getADeriv(self, u, v, adjoint= False):
+
+        D = self.Div
+        G = self.Grad
+        MfRhoIDeriv = self.MfRhoIDeriv
+
+        if adjoint:
+            # if self._makeASymmetric is True:
+            #     v = V * v
+            return(MfRhoIDeriv( G * u ).T) * ( D.T * v)
+
+        # I think we should deprecate this for DC problem.
+        # if self._makeASymmetric is True:
+        #     return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) )
+        return D * (MfRhoIDeriv( G * u ) * v)
+
+    def getRHS(self):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm()
+
+        # I think we should deprecate this for DC problem.
+        # if self._makeASymmetric is True:
+        #     return self.Vol.T * RHS
+
+        return RHS
+
+    def getRHSDeriv(self, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
+        # return qDeriv
+        return Zero()
+
+    def setBC(self):
+        if self.mesh.dim==3:
+            fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+            gBFzm = self.mesh.gridFz[fzm,:]
+            gBFzp = self.mesh.gridFz[fzp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+            temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+            alpha_zm, alpha_zp = temp_zm*0., temp_zp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+            beta_zm, beta_zp = temp_zm, temp_zp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+            gamma_zm, gamma_zp = temp_zm*0., temp_zp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
+
+        elif self.mesh.dim==2:
+
+            fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
+
+        x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma)
+        V = self.Vol
+        self.Div = V * self.mesh.faceDiv
+        P_BC, B = self.mesh.getBCProjWF_simple()
+        M = B*self.mesh.aveCC2F
+        self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
+
+
+class Problem3D_N(BaseIPProblem):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'EB' # N potentials means B is on faces
+    fieldsPair    = Fields_N
+
+    def __init__(self, mesh, **kwargs):
+        BaseIPProblem.__init__(self, mesh, **kwargs)
+
+    def getA(self):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = G.T MeSigma G
+
+        """
+
+        MeSigma = self.MeSigma
+        Grad = self.mesh.nodalGrad
+        A = Grad.T * MeSigma * Grad
+
+        # Handling Null space of A
+        A[0,0] = A[0,0] + 1.
+
+        return A
+
+    def getADeriv(self, u, v, adjoint=False):
+        """
+
+        Product of the derivative of our system matrix with respect to the model and a vector
+
+        """
+        MeSigma = self.MeSigma
+        Grad = self.mesh.nodalGrad
+        if not adjoint:
+            return Grad.T*(self.MeSigmaDeriv(Grad*u)*v)
+        elif adjoint:
+            return self.MeSigmaDeriv(Grad*u).T * (Grad*v)
+
+
+    def getRHS(self):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm()
+        return RHS
+
+    def getRHSDeriv(self, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
+        # return qDeriv
+        return Zero()
+
+if __name__ == '__main__':
+
+
+    cs = 12.5
+    hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
+    hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
+    hz = [(cs,7, -1.3),(cs,20)]
+    mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
+    sigma = np.ones(mesh.nC)
+    prob = BaseIPProblem(mesh, sigma=sigma)
+
+

SimPEG/EM/Static/IP/SurveyIP.py

@@ -0,0 +1,23 @@
+import SimPEG
+from SimPEG.EM.Base import BaseEMSurvey
+from SimPEG import sp, Survey
+from SimPEG.Utils import Zero, Identity
+from SimPEG.EM.Static.DC.SrcDC import BaseSrc
+from SimPEG.EM.Static.DC.RxDC import BaseRx
+
+class Survey(BaseEMSurvey):
+    rxPair  = BaseRx
+    srcPair = BaseSrc
+
+    def __init__(self, srcList, **kwargs):
+        self.srcList = srcList
+        BaseEMSurvey.__init__(self, srcList, **kwargs)
+
+    def dpred(self, m, f=None):
+        """
+            Predicted data.
+
+            .. math::
+                d_\\text{pred} = Pf(m)
+        """
+        return self.prob.Jvec(m, m, f=f)

SimPEG/EM/Static/IP/__init__.py

@@ -0,0 +1,2 @@
+from ProblemIP import Problem3D_CC, Problem3D_N
+from SurveyIP import Survey

SimPEG/EM/Static/SIP/ProblemSIP.py

@@ -0,0 +1,445 @@
+from SimPEG import Problem, Utils, Maps, Mesh
+from SimPEG.EM.Base import BaseEMProblem
+from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
+from SimPEG.Utils import sdiag
+import numpy as np
+from SimPEG.Utils import Zero
+from SimPEG.EM.Static.DC import getxBCyBC_CC
+from SurveySIP import Survey, Data
+
+class ColeColePropMap(Maps.PropMap):
+    """
+        Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m)
+    """
+
+    eta = Maps.Property("Electrical Conductivity", defaultInvProp=True)
+    tau = Maps.Property("Electrical Conductivity", defaultVal=0.1, propertyLink=('taui', Maps.ReciprocalMap))
+    taui   = Maps.Property("Electrical Conductivity", defaultVal=1., propertyLink=('tau', Maps.ReciprocalMap))
+    c   = Maps.Property("Electrical Conductivity", defaultVal=1.)
+
+
+class BaseSIPProblem(BaseEMProblem):
+
+    surveyPair = Survey
+    fieldsPair = Fields
+    dataPair = Data
+    PropMap = ColeColePropMap
+    Ainv = None
+    sigma = None
+    rho = None
+    f = None
+    Ainv = None
+
+    def DebyeTime(self, t):
+        peta = self.curModel.eta*np.exp(-self.curModel.taui*t)
+        return peta
+
+    def EtaDeriv(self, t, v, adjoint=False):
+        v = np.array(v, dtype=float)
+        if adjoint:
+            return self.curModel.etaDeriv.T * (np.exp(-self.curModel.taui*t)*v)
+        else:
+            return np.exp(-self.curModel.taui*t) * (self.curModel.etaDeriv*v)
+
+
+    def TauiDeriv(self, t, v, adjoint=False):
+        v = np.array(v, dtype=float)
+        if adjoint:
+            return -self.curModel.tauiDeriv.T * (self.curModel.eta*t*np.exp(-self.curModel.taui*t)*v)
+        else:
+            return -self.curModel.eta*t*np.exp(-self.curModel.taui*t) * (self.curModel.tauiDeriv*v)
+
+    def fields(self, m):
+        self.curModel = m
+        if self.f is None:
+            self.f = self.fieldsPair(self.mesh, self.survey)
+            if self.Ainv == None:
+                A = self.getA()
+                self.Ainv = self.Solver(A, **self.solverOpts)
+            RHS = self.getRHS()
+            u = self.Ainv * RHS
+            Srcs = self.survey.srcList
+            self.f[Srcs, self._solutionType] = u
+        return self.f
+
+    def forward(self, m, f=None):
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+        Jv = self.dataPair(self.survey) #same size as the data
+        # A = self.getA()
+        JvAll = []
+        for tind in range(len(self.survey.times)):
+            #Pseudo-chareability
+            t = self.survey.times[tind]
+            v = self.DebyeTime(t)
+            for src in self.survey.srcList:
+                u_src = f[src, self._solutionType] # solution vector
+                dA_dm_v = self.getADeriv(u_src, v)
+                dRHS_dm_v = self.getRHSDeriv(src, v)
+                du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
+                for rx in src.rxList:
+                    timeindex = rx.getTimeP(self.survey.times)
+                    if timeindex[tind]:
+                        df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                        df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
+                        Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
+
+        # Conductivity (d u / d log sigma)
+        if self._formulation is 'EB':
+            return -Utils.mkvc(Jv)
+        # Resistivity (d u / d log rho)
+        if self._formulation is 'HJ':
+            return Utils.mkvc(Jv)
+
+    def Jvec(self, m, v, f=None):
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+        Jv = self.dataPair(self.survey) #same size as the data
+        # A = self.getA()
+        JvAll = []
+        #Assume only eta and tau (eta first then tau)
+        # v = [2*Mx1]
+        v = v.reshape((int(v.size/2), 2), order='F')
+
+        for tind in range(len(self.survey.times)):
+            t = self.survey.times[tind]
+            v0 = self.EtaDeriv(t, v[:,0])
+            v1 = self.TauiDeriv(t, v[:,1])
+            for src in self.survey.srcList:
+                u_src = f[src, self._solutionType] # solution vector
+                dA_dm_v0 = self.getADeriv(u_src, v0)
+                dRHS_dm_v0 = self.getRHSDeriv(src, v0)
+                du_dm_v0 = self.Ainv * ( - dA_dm_v0 + dRHS_dm_v0 )
+                dA_dm_v1 = self.getADeriv(u_src, v1)
+                dRHS_dm_v1 = self.getRHSDeriv(src, v1)
+                du_dm_v1 = self.Ainv * ( - dA_dm_v1 + dRHS_dm_v1 )
+                for rx in src.rxList:
+                    timeindex = rx.getTimeP(self.survey.times)
+                    if timeindex[tind]:
+                        df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                        df_dm_v0 = df_dmFun(src, du_dm_v0, v0, adjoint=False)
+                        df_dm_v1 = df_dmFun(src, du_dm_v1, v1, adjoint=False)
+                        Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v0)
+                        Jv[src, rx, t] += rx.evalDeriv(src, self.mesh, f, df_dm_v1)
+        # Conductivity (d u / d log sigma)
+        if self._formulation is 'EB':
+            return -Jv.tovec()
+        # Resistivity (d u / d log rho)
+        if self._formulation is 'HJ':
+            return Jv.tovec()
+
+    def Jtvec(self, m, v, f=None):
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        # Ensure v is a data object.
+        if not isinstance(v, self.dataPair):
+            v = self.dataPair(self.survey, v)
+
+        Jtv= np.zeros(m.size)
+        for tind in range(len(self.survey.times)):
+            t = self.survey.times[tind]
+            for src in self.survey.srcList:
+                u_src = f[src, self._solutionType]
+                for rx in src.rxList:
+                    timeindex = rx.getTimeP(self.survey.times)
+                    if timeindex[tind]:
+                        PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx, t], adjoint=True) # wrt f, need possibility wrt m
+                        df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                        df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
+                        ATinvdf_duT = self.Ainv * df_duT
+                        dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
+                        dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
+                        du_dmT = -dA_dmT + dRHS_dmT
+                        Jtv += np.r_[self.EtaDeriv(self.survey.times[tind], du_dmT, adjoint=True), self.TauiDeriv(self.survey.times[tind], du_dmT, adjoint=True)]
+
+        # Conductivity ((d u / d log sigma).T)
+        if self._formulation is 'EB':
+            return -Jtv
+        # Conductivity ((d u / d log rho).T)
+        if self._formulation is 'HJ':
+            return Jtv
+
+    def getSourceTerm(self):
+        """
+        takes concept of source and turns it into a matrix
+        """
+        """
+        Evaluates the sources, and puts them in matrix form
+
+        :rtype: (numpy.ndarray, numpy.ndarray)
+        :return: q (nC or nN, nSrc)
+        """
+
+        Srcs = self.survey.srcList
+
+        if self._formulation is 'EB':
+            n = self.mesh.nN
+            # return NotImplementedError
+
+        elif self._formulation is 'HJ':
+            n = self.mesh.nC
+
+        q = np.zeros((n, len(Srcs)))
+
+        for i, src in enumerate(Srcs):
+            q[:,i] = src.eval(self)
+        return q
+
+    @property
+    def deleteTheseOnModelUpdate(self):
+        toDelete = []
+        return toDelete
+
+    # assume log rho or log cond
+    @property
+    def MeSigma(self):
+        """
+            Edge inner product matrix for \\(\\sigma\\). Used in the E-B formulation
+        """
+        if getattr(self, '_MeSigma', None) is None:
+            self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma)
+        return self._MeSigma
+
+    @property
+    def MfRhoI(self):
+        """
+            Inverse of :code:`MfRho`
+        """
+        if getattr(self, '_MfRhoI', None) is None:
+            self._MfRhoI = self.mesh.getFaceInnerProduct(self.rho, invMat=True)
+        return self._MfRhoI
+
+    def MfRhoIDeriv(self,u):
+        """
+            Derivative of :code:`MfRhoI` with respect to the model.
+        """
+
+        dMfRhoI_dI = -self.MfRhoI**2
+        dMf_drho = self.mesh.getFaceInnerProductDeriv(self.rho)(u)
+        drho_dlogrho = Utils.sdiag(self.rho)
+        return dMfRhoI_dI * ( dMf_drho * ( drho_dlogrho))
+
+    # TODO: This should take a vector
+    def MeSigmaDeriv(self, u):
+        """
+            Derivative of MeSigma with respect to the model
+        """
+        dsigma_dlogsigma = Utils.sdiag(self.sigma)
+        return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) * dsigma_dlogsigma
+
+class Problem3D_CC(BaseSIPProblem):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'HJ' # CC potentials means J is on faces
+    fieldsPair    = Fields_CC
+
+    def __init__(self, mesh, **kwargs):
+        BaseSIPProblem.__init__(self, mesh, **kwargs)
+        self.setBC()
+
+    def getA(self):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = D MfRhoI G
+
+        """
+
+        D = self.Div
+        G = self.Grad
+        # TODO: this won't work for full anisotropy
+        MfRhoI = self.MfRhoI
+        A = D * MfRhoI * G
+
+        # I think we should deprecate this for DC problem.
+        # if self._makeASymmetric is True:
+        #     return V.T * A
+        return A
+
+    def getADeriv(self, u, v, adjoint= False):
+
+        D = self.Div
+        G = self.Grad
+        MfRhoIDeriv = self.MfRhoIDeriv
+
+        if adjoint:
+            # if self._makeASymmetric is True:
+            #     v = V * v
+            return(MfRhoIDeriv( G * u ).T) * ( D.T * v)
+
+        # I think we should deprecate this for DC problem.
+        # if self._makeASymmetric is True:
+        #     return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) )
+        return D * (MfRhoIDeriv( G * u ) * v)
+
+    def getRHS(self):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm()
+
+        # I think we should deprecate this for DC problem.
+        # if self._makeASymmetric is True:
+        #     return self.Vol.T * RHS
+
+        return RHS
+
+    def getRHSDeriv(self, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
+        # return qDeriv
+        return Zero()
+
+    def setBC(self):
+        if self.mesh.dim==3:
+            fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+            gBFzm = self.mesh.gridFz[fzm,:]
+            gBFzp = self.mesh.gridFz[fzp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+            temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+            alpha_zm, alpha_zp = temp_zm*0., temp_zp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+            beta_zm, beta_zp = temp_zm, temp_zp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+            gamma_zm, gamma_zp = temp_zm*0., temp_zp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
+
+        elif self.mesh.dim==2:
+
+            fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
+
+        x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma)
+        V = self.Vol
+        self.Div = V * self.mesh.faceDiv
+        P_BC, B = self.mesh.getBCProjWF_simple()
+        M = B*self.mesh.aveCC2F
+        self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
+
+
+class Problem3D_N(BaseSIPProblem):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'EB' # N potentials means B is on faces
+    fieldsPair    = Fields_N
+
+    def __init__(self, mesh, **kwargs):
+        BaseSIPProblem.__init__(self, mesh, **kwargs)
+
+    def getA(self):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = G.T MeSigma G
+
+        """
+
+        # TODO: this won't work for full anisotropy
+        MeSigma = self.MeSigma
+        Grad = self.mesh.nodalGrad
+        A = Grad.T * MeSigma * Grad
+
+        # Handling Null space of A
+        A[0,0] = A[0,0] + 1.
+
+        return A
+
+    def getADeriv(self, u, v, adjoint=False):
+        """
+
+        Product of the derivative of our system matrix with respect to the model and a vector
+
+        """
+        MeSigma = self.MeSigma
+        Grad = self.mesh.nodalGrad
+        if not adjoint:
+            return Grad.T*(self.MeSigmaDeriv(Grad*u)*v)
+        elif adjoint:
+            return self.MeSigmaDeriv(Grad*u).T * (Grad*v)
+
+
+    def getRHS(self):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm()
+        return RHS
+
+    def getRHSDeriv(self, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
+        # return qDeriv
+        return Zero()
+
+if __name__ == '__main__':
+
+
+    cs = 12.5
+    hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
+    hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
+    hz = [(cs,7, -1.3),(cs,20)]
+    mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
+    sigma = np.ones(mesh.nC)
+    prob = BaseSIPProblem(mesh, sigma=sigma)
+
+

SimPEG/EM/Static/SIP/Regularization.py

@@ -0,0 +1,204 @@
+from SimPEG import Utils, Maps, Mesh, sp, np
+from SimPEG.Regularization import BaseRegularization, Simple
+
+class MultiRegularization(Simple):
+    """
+    **MultiRegularization Class**
+
+    This is used to regularize the model space
+    having multiple models [m1, m2, m3, ...] ::
+
+        reg = Regularization(mesh)
+
+    """
+    nModels = None # Number of models
+    ratios = None  # Ratio for different models
+    crossgrad = False # Use cross gradient or not
+    betacross = 1.
+    wx = []
+    wy = []
+    wz = []
+
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        if self.nModels == None:
+            raise Exception("Put nModels as a initial input!")
+        if self.ratios == None:
+            self.ratios = [1. for imodel in range(self.nModels)]
+
+    @property
+    def Wsmall(self):
+        """Regularization matrix Wsmall"""
+        if getattr(self,'_Wsmall', None) is None:
+            vecs = []
+            for imodel in range(self.nModels):
+                vecs.append((self.regmesh.vol*self.alpha_s*self.wght*self.ratios[imodel])**0.5)
+            self._Wsmall = Utils.sdiag(np.hstack(vecs))
+        return self._Wsmall
+
+    @property
+    def Wx(self):
+        """Regularization matrix Wx"""
+        if getattr(self, '_Wx', None) is None:
+            mats = []
+            for imodel in range(self.nModels):
+                self.wx.append(Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*self.ratios[imodel]*(self.regmesh.aveCC2Fx*self.wght))**0.5))
+                mats.append(self.wx[imodel]*self.regmesh.cellDiffxStencil)
+            self._Wx = sp.block_diag(mats)
+        return self._Wx
+
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
+        if getattr(self, '_Wy', None) is None:
+            mats = []
+            for imodel in range(self.nModels):
+                self.wy.append(Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol*self.alpha_y*self.ratios[imodel]*(self.regmesh.aveCC2Fy*self.wght))**0.5))
+                mats.append(self.wy[imodel]*self.regmesh.cellDiffyStencil)
+            self._Wy = sp.block_diag(mats)
+        return self._Wy
+
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
+        if getattr(self, '_Wz', None) is None:
+            mats = []
+            for imodel in range(self.nModels):
+                self.wz.append(Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*self.ratios[imodel]*(self.regmesh.aveCC2Fz*self.wght))**0.5))
+                mats.append(self.wz[imodel]*self.regmesh.cellDiffzStencil)
+            self._Wz = sp.block_diag(mats)
+        return self._Wz
+
+    @property
+    def Wsmooth(self):
+        """Full smoothness regularization matrix W"""
+        if getattr(self, '_Wsmooth', None) is None:
+            wlist = (self.Wx,)
+            if self.regmesh.dim > 1:
+                wlist += (self.Wy,)
+            if self.regmesh.dim > 2:
+                wlist += (self.Wz,)
+            self._Wsmooth = sp.vstack(wlist)
+        return self._Wsmooth
+
+    @property
+    def W(self):
+        """Full regularization matrix W"""
+        if getattr(self, '_W', None) is None:
+            wlist = (self.Wsmall, self.Wsmooth)
+            self._W = sp.vstack(wlist)
+        return self._W
+
+
+    @Utils.timeIt
+    def eval(self, m):
+        return self._evalSmall(m) + self._evalSmooth(m)
+
+    @Utils.timeIt
+    def _evalSmall(self, m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmooth(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wsmooth * ( self.mapping * m)
+        return 0.5 * r.dot(r)
+
+    def cross(a,b):
+        ax, ay, az = a[0], a[1], a[2]
+        bx, by, bz = b[0], b[1], b[2]
+        cx = ay*bz - az*by
+        cy = az*bx - ax*bz
+        cz = ax*by - ay*bx
+        return [cx, cy, cz]
+
+     # TODO: Implement Cross Gradients..
+    @Utils.timeIt
+    def _evalCross(self, m):
+        if self.crossgrad == False:
+            return 0.
+        elif self.crossgrad == True:
+            M = (self.mapping * m).reshape((self.regmesh.nC, self.nModels), order="F")
+
+            ax = self.regmesh.aveFx2CC*self.regmesh.wx[0]*M[:,0]
+            ay = self.regmesh.aveFy2CC*self.regmesh.wy[0]*M[:,0]
+            az = self.regmesh.aveFz2CC*self.regmesh.wz[0]*M[:,0]
+            bx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1]
+            by = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1]
+            bz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1]
+            #ab
+            out_ab = cross([ax, ay, az], [bx, by, bz])
+            r = np.r_[out_ab[0], out_ab[1], out_ab[2]]*np.sqrt(self.betacross)
+
+            if self.nModels == 3:
+                cx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1]
+                cy = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1]
+                cz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1]
+                #ac
+                out_ac = cross([ax, ay, az], [cx, cy, cz])
+                #bc
+                out_bc = cross([bx, by, bz], [cx, cy, cz])
+                r = np.r_[r, np.hstack(out_ac)*np.sqrt(self.betacross), np.hstack(out_bc)*np.sqrt(self.betacross)]
+
+            return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def evalDeriv(self, m):
+        """
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        deriv = self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+        if self.crossgrad==True:
+            deriv += self._evalCrossDeriv(m)
+        return deriv
+
+    @Utils.timeIt
+    def _evalCrossDeriv(self,m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
+        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
+
+    @Utils.timeIt
+    def eval2Deriv(self, m, v=None):
+        """
+        Second derivative
+
+        :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::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the second derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W}
+
+        """
+        mD = self.mapping.deriv(m - self.mref)
+        if v is None:
+            return mD.T * self.W.T * self.W * mD
+
+        return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
+
+
+

SimPEG/EM/Static/SIP/RxSIP.py

@@ -0,0 +1,88 @@
+import SimPEG
+import numpy as np
+from SimPEG.Utils import Zero, closestPoints
+
+class BaseRx(SimPEG.Survey.BaseTimeRx):
+    locs = None
+    rxType = None
+
+    knownRxTypes = {
+                    'phi':['phi',None],
+                    'ex':['e','x'],
+                    'ey':['e','y'],
+                    'ez':['e','z'],
+                    'jx':['j','x'],
+                    'jy':['j','y'],
+                    'jz':['j','z'],
+                    }
+
+    def __init__(self, locs, times, rxType, **kwargs):
+        SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType, **kwargs)
+
+    @property
+    def projField(self):
+        """Field Type projection (e.g. e b ...)"""
+        return self.knownRxTypes[self.rxType][0]
+
+    def projGLoc(self, f):
+        """Grid Location projection (e.g. Ex Fy ...)"""
+        comp = self.knownRxTypes[self.rxType][1]
+        if comp is not None:
+            return f._GLoc(self.rxType) + comp
+        return f._GLoc(self.rxType)
+
+    def getTimeP(self, timesall):
+        """
+            Returns the time projection matrix.
+
+            .. note::
+
+                This is not stored in memory, but is created on demand.
+        """
+        time_inds = np.in1d(timesall, self.times)
+        return time_inds
+
+    def evalDeriv(self, src, mesh, f, v, adjoint=False):
+        P = self.getP(mesh, self.projGLoc(f))
+        if not adjoint:
+            return P*v
+        elif adjoint:
+            return P.T*v
+
+
+# DC.Rx.Dipole(locs)
+class Dipole(BaseRx):
+
+    def __init__(self, locsM, locsN, times, rxType = 'phi', **kwargs):
+        assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size'
+        locs = [locsM, locsN]
+        # We may not need this ...
+        BaseRx.__init__(self, locs, times, rxType)
+
+    @property
+    def nD(self):
+        """Number of data in the receiver."""
+        # return self.locs[0].shape[0] * len(self.times)
+        return self.locs[0].shape[0]
+
+    @property
+    def nRx(self):
+        """Number of data in the receiver."""
+        return self.locs[0].shape[0]
+
+        # Not sure why ...
+        # return int(self.locs[0].size / 2)
+
+
+    def getP(self, mesh, Gloc):
+        if mesh in self._Ps:
+            return self._Ps[mesh]
+
+        P0 = mesh.getInterpolationMat(self.locs[0], Gloc)
+        P1 = mesh.getInterpolationMat(self.locs[1], Gloc)
+        P = P0 - P1
+
+        if self.storeProjections:
+            self._Ps[mesh] = P
+
+        return P

SimPEG/EM/Static/SIP/SrcSIP.py

@@ -0,0 +1,64 @@
+import SimPEG
+# from SimPEG.EM.Base import BaseEMSurvey
+from SimPEG.Utils import Zero, closestPoints, mkvc
+import numpy as np
+
+class BaseSrc(SimPEG.Survey.BaseSrc):
+
+    current = 1.0
+    loc = None
+
+    def __init__(self, rxList, **kwargs):
+        SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs)
+
+    def eval(self, prob):
+        raise NotImplementedError
+
+    def evalDeriv(self, prob):
+        return Zero()
+
+    @property
+    def nD(self):
+        """Number of data"""
+        return self.vnD.sum()
+
+    @property
+    def vnD(self):
+        """Vector number of data"""
+        return np.array([rx.nD*len(rx.times) for rx in self.rxList])
+
+
+
+class Dipole(BaseSrc):
+
+    def __init__(self, rxList, locA, locB, **kwargs):
+        assert locA.shape == locB.shape, 'Shape of locA and locB should be the same'
+        self.loc = [locA, locB]
+        BaseSrc.__init__(self, rxList, **kwargs)
+
+    def eval(self, prob):
+        if prob._formulation == 'HJ':
+            inds = closestPoints(prob.mesh, self.loc, gridLoc='CC')
+            q = np.zeros(prob.mesh.nC)
+            q[inds] = self.current * np.r_[1., -1.]
+        elif prob._formulation == 'EB':
+            qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense()
+            qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense()
+            q = self.current * mkvc(qa+qb)
+        return q
+
+class Pole(BaseSrc):
+
+    def __init__(self, rxList, loc, **kwargs):
+        BaseSrc.__init__(self, rxList, loc=loc, **kwargs)
+
+    def eval(self, prob):
+        if prob._formulation == 'HJ':
+            inds = closestPoints(prob.mesh, self.loc)
+            q = np.zeros(prob.mesh.nC)
+            q[inds] = self.current * np.r_[1.]
+        elif prob._formulation == 'EB':
+            q = prob.mesh.getInterpolationMat(self.loc, locType='N').todense()
+            q = self.current * mkvc(q)
+        return q
+

SimPEG/EM/Static/SIP/SurveySIP.py

@@ -0,0 +1,102 @@
+import SimPEG
+from SimPEG.EM.Base import BaseEMSurvey
+from SimPEG import np, sp, Survey, Utils
+from SimPEG.Utils import Zero, Identity
+from SimPEG.EM.Static.SIP.SrcSIP import BaseSrc
+from SimPEG.EM.Static.SIP.RxSIP import BaseRx
+import uuid
+
+
+class Survey(BaseEMSurvey):
+    rxPair  = BaseRx
+    srcPair = BaseSrc
+    times = None
+
+    def __init__(self, srcList, **kwargs):
+        self.srcList = srcList
+        BaseEMSurvey.__init__(self, srcList, **kwargs)
+        self.getUniqueTimes()
+
+    def getUniqueTimes(self):
+        time_rx = []
+        for src in self.srcList:
+            for rx in src.rxList:
+                time_rx.append(rx.times)
+        self.times = np.unique(np.hstack(time_rx))
+
+    def dpred(self, m, f=None):
+        """
+            Predicted data.
+
+            .. math::
+                d_\\text{pred} = Pf(m)
+        """
+        return self.prob.forward(m, f=f)
+
+
+class Data(SimPEG.Survey.Data):
+    """Fancy data storage by Src and Rx"""
+
+    def __init__(self, survey, v=None):
+        self.uid = str(uuid.uuid4())
+        self.survey = survey
+        self._dataDict = {}
+        for src in self.survey.srcList:
+            self._dataDict[src] = {}
+            for rx in src.rxList:
+                self._dataDict[src][rx] = {}
+
+        if v is not None:
+            self.fromvec(v)
+
+    def _ensureCorrectKey(self, key):
+        if type(key) is tuple:
+            if len(key) is not 3:
+                raise KeyError('Key must be [Src, Rx, tInd]')
+            if key[0] not in self.survey.srcList:
+                raise KeyError('Src Key must be a source in the survey.')
+            if key[1] not in key[0].rxList:
+                raise KeyError('Rx Key must be a receiver for the source.')
+            return key
+        elif isinstance(key, self.survey.srcPair):
+            if key not in self.survey.srcList:
+                raise KeyError('Key must be a source in the survey.')
+            return key, None, None
+        else:
+            raise KeyError('Key must be [Src] or [Src,Rx] or [Src, Rx, tInd]')
+
+    def __setitem__(self, key, value):
+        src, rx, t = self._ensureCorrectKey(key)
+        assert rx is not None, 'set data using [Src, Rx]'
+        assert isinstance(value, np.ndarray), 'value must by ndarray'
+        assert value.size == rx.nD, "value must have the same number of data as the source."
+        self._dataDict[src][rx][t] = Utils.mkvc(value)
+
+    def __getitem__(self, key):
+        src, rx, t = self._ensureCorrectKey(key)
+        if rx is not None:
+            if rx not in self._dataDict[src]:
+                raise Exception('Data for receiver has not yet been set.')
+            return self._dataDict[src][rx][t]
+
+        return np.concatenate([self[src,rx, t] for rx in src.rxList])
+
+    def tovec(self):
+        val = []
+        for src in self.survey.srcList:
+            for rx in src.rxList:
+                for t in rx.times:
+                    val.append(self[src, rx, t])
+        return np.concatenate(val)
+
+
+    def fromvec(self, v):
+        v = Utils.mkvc(v)
+        assert v.size == self.survey.nD, 'v must have the correct number of data.'
+        indBot, indTop = 0, 0
+        for src in self.survey.srcList:
+            for rx in src.rxList:
+                for t in rx.times:
+                    indTop += rx.nRx
+                    self[src, rx, t] = v[indBot:indTop]
+                    indBot += rx.nRx

SimPEG/EM/Static/SIP/__init__.py

@@ -0,0 +1,5 @@
+from ProblemSIP import Problem3D_CC, Problem3D_N
+from SurveySIP import Survey, Data
+import SrcSIP as Src #Pole
+import RxSIP as Rx
+from Regularization import MultiRegularization

SimPEG/EM/Static/Utils/StaticUtils.py

@@ -0,0 +1,421 @@
+from SimPEG import np
+from SimPEG.EM.Static import DC, IP
+
+def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
+    """
+        Read list of 2D tx-rx location and plot a speudo-section of apparent
+        resistivity.
+
+        Assumes flat topo for now...
+
+        Input:
+        :param d2D, z0
+        :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole)
+        :switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential)
+        Output:
+        :figure scatter plot overlayed on image
+
+        Edited Feb 17th, 2016
+
+        @author: dominiquef
+
+    """
+    from SimPEG import np
+    from scipy.interpolate import griddata
+    import pylab as plt
+
+    # Set depth to 0 for now
+    z0 = 0.
+
+    # Pre-allocate
+    midx = []
+    midz = []
+    rho = []
+    LEG = []
+    count = 0 # Counter for data
+    for ii in range(DCsurvey.nSrc):
+
+        Tx = DCsurvey.srcList[ii].loc
+        Rx = DCsurvey.srcList[ii].rxList[0].locs
+
+        nD = DCsurvey.srcList[ii].rxList[0].nD
+
+        data = DCsurvey.dobs[count:count+nD]
+        count += nD
+
+        # Get distances between each poles A-B-M-N
+        if stype == 'pdp':
+            MA = np.abs(Tx[0] - Rx[0][:,0])
+            NA = np.abs(Tx[0] - Rx[1][:,0])
+            MN = np.abs(Rx[1][:,0] - Rx[0][:,0])
+
+            # Create mid-point location
+            Cmid = Tx[0]
+            Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
+            if DCsurvey.mesh.dim == 2:
+                zsrc = Tx[1]
+            elif DCsurvey.mesh.dim ==3:
+                zsrc = Tx[2]
+
+        elif stype == 'dpdp':
+            MA = np.abs(Tx[0][0] - Rx[0][:,0])
+            MB = np.abs(Tx[1][0] - Rx[0][:,0])
+            NA = np.abs(Tx[0][0] - Rx[1][:,0])
+            NB = np.abs(Tx[1][0] - Rx[1][:,0])
+
+            # Create mid-point location
+            Cmid = (Tx[0][0] + Tx[1][0])/2
+            Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
+            if DCsurvey.mesh.dim == 2:
+                zsrc = (Tx[0][1] + Tx[1][1])/2
+            elif DCsurvey.mesh.dim ==3:
+                zsrc = (Tx[0][2] + Tx[1][2])/2
+
+        # Change output for dtype
+        if dtype == 'volt':
+
+            rho = np.hstack([rho,data])
+
+        else:
+
+            # Compute pant leg of apparent rho
+            if stype == 'pdp':
+
+                leg =  data * 2*np.pi  * MA * ( MA + MN ) / MN
+
+            elif stype == 'dpdp':
+
+                leg = data * 2*np.pi / ( 1/MA - 1/MB + 1/NB - 1/NA )
+                LEG.append(1./(2*np.pi) *( 1/MA - 1/MB + 1/NB - 1/NA ))
+            else:
+                print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
+                break
+
+
+            if dtype == 'appc':
+
+                leg = np.log10(abs(1./leg))
+                rho = np.hstack([rho,leg])
+
+            elif dtype == 'appr':
+
+                leg = np.log10(abs(leg))
+                rho = np.hstack([rho,leg])
+
+            else:
+                print """dtype must be 'appr' | 'appc' | 'volt' """
+                break
+
+
+        midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
+        if DCsurvey.mesh.dim==3:
+            midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ])
+        elif DCsurvey.mesh.dim==2:
+            midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ])
+    ax = axs
+
+    # Grid points
+    grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
+    grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear')
+
+    if clim == None:
+        vmin, vmax = rho.min(), rho.max()
+    else:
+        vmin, vmax = clim[0], clim[1]
+
+    grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho)
+    ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, clim=(vmin, vmax), vmin=vmin, vmax=vmax)
+    cbar = plt.colorbar(format="$10^{%.1f}$",fraction=0.04,orientation="horizontal")
+
+    cmin,cmax = cbar.get_clim()
+    ticks = np.linspace(cmin,cmax,3)
+    cbar.set_ticks(ticks)
+    cbar.ax.tick_params(labelsize=10)
+
+    if dtype == 'appc':
+        cbar.set_label("App.Cond",size=12)
+    elif dtype == 'appr':
+        cbar.set_label("App.Res.",size=12)
+    elif dtype == 'volt':
+        cbar.set_label("Potential (V)",size=12)
+
+    # Plot apparent resistivity
+    ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
+
+    #ax.set_xticklabels([])
+    #ax.set_yticklabels([])
+
+    plt.gca().set_aspect('equal', adjustable='box')
+
+
+
+    return ph, LEG
+
+def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
+    """
+        Load in endpoints and survey specifications to generate Tx, Rx location
+        stations.
+
+        Assumes flat topo for now...
+
+        Input:
+        :param endl -> input endpoints [x1, y1, z1, x2, y2, z2]
+        :object mesh -> SimPEG mesh object
+        :switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient'
+        : param a, n -> pole seperation, number of rx dipoles per tx
+
+        Output:
+        :param Tx, Rx -> List objects for each tx location
+            Lines: P1x, P1y, P1z, P2x, P2y, P2z
+
+        Created on Wed December 9th, 2015
+
+        @author: dominiquef
+        !! Require clean up to deal with DCsurvey
+    """
+
+    from SimPEG import np
+
+    def xy_2_r(x1,x2,y1,y2):
+        r = np.sqrt( np.sum((x2 - x1)**2 + (y2 - y1)**2) )
+        return r
+
+    ## Evenly distribute electrodes and put on surface
+    # Mesure survey length and direction
+    dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
+
+    dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
+    dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
+
+    nstn = np.floor( dl_len / a )
+
+    # Compute discrete pole location along line
+    stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
+    stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
+
+    if mesh.dim==2:
+        ztop = mesh.vectorNy[-1]
+        # Create line of P1 locations
+        M = np.c_[stn_x, np.ones(nstn).T*ztop]
+        # Create line of P2 locations
+        N = np.c_[stn_x+a*dl_x, np.ones(nstn).T*ztop]
+
+    elif mesh.dim==3:
+        ztop = mesh.vectorNz[-1]
+        # Create line of P1 locations
+        M = np.c_[stn_x, stn_y, np.ones(nstn).T*ztop]
+        # Create line of P2 locations
+        N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*ztop]
+
+
+    ## Build list of Tx-Rx locations depending on survey type
+    # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn]
+    # Pole-dipole: Moving pole on one end -> [A a MN1 a MN2 ... MNn a B]
+    SrcList = []
+
+
+    if stype != 'gradient':
+
+        for ii in range(0, int(nstn)-1):
+
+
+            if stype == 'dpdp':
+                tx = np.c_[M[ii,:],N[ii,:]]
+            elif stype == 'pdp':
+                tx = np.c_[M[ii,:],M[ii,:]]
+
+            # Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]])
+
+            # Current elctrode seperation
+            AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
+
+            # Number of receivers to fit
+            nstn = np.min([np.floor( (AB - b) / a ) , n])
+
+            # Check if there is enough space, else break the loop
+            if nstn <= 0:
+                continue
+
+            # Compute discrete pole location along line
+            stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
+            stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
+
+            # Create receiver poles
+
+            if mesh.dim==3:
+                # Create line of P1 locations
+                P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*ztop]
+                # Create line of P2 locations
+                P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*ztop]
+                rxClass = DC.Rx.Dipole(P1, P2)
+
+            elif mesh.dim==2:
+                # Create line of P1 locations
+                P1 = np.c_[stn_x, np.ones(nstn).T*ztop]
+                # Create line of P2 locations
+                P2 = np.c_[stn_x+a*dl_x, np.ones(nstn).T*ztop]
+                rxClass = DC.Rx.Dipole_ky(P1, P2)
+
+            if stype == 'dpdp':
+                    srcClass = DC.Src.Dipole([rxClass], M[ii,:],N[ii,:])
+            elif stype == 'pdp':
+                    srcClass = DC.Src.Pole([rxClass], M[ii,:])
+            SrcList.append(srcClass)
+
+    elif stype == 'gradient':
+
+        # Gradient survey only requires Tx at end of line and creates a square
+        # grid of receivers at in the middle at a pre-set minimum distance
+
+        # Get the edge limit of survey area
+        min_x = endl[0,0] + dl_x * b
+        min_y = endl[0,1] + dl_y * b
+
+        max_x = endl[1,0] - dl_x * b
+        max_y = endl[1,1] - dl_y * b
+
+        box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
+        box_w = box_l/2.
+
+        nstn = np.floor( box_l / a )
+
+        # Compute discrete pole location along line
+        stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
+        stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
+
+        # Define number of cross lines
+        nlin = int(np.floor( box_w / a ))
+        lind = range(-nlin,nlin+1)
+
+        ngrad = nstn * len(lind)
+
+        rx = np.zeros([ngrad,6])
+        for ii in range( len(lind) ):
+
+            # Move line in perpendicular direction by dipole spacing
+            lxx = stn_x - lind[ii]*a*dl_y
+            lyy = stn_y + lind[ii]*a*dl_x
+
+
+            M = np.c_[ lxx, lyy , np.ones(nstn).T*ztop]
+            N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*ztop]
+            rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N]
+
+            if mesh.dim==3:
+                rxClass = DC.Rx.Dipole(rx[:,:3], rx[:,3:])
+            elif mesh.dim==2:
+                M = M[:,[0,2]]
+                N = N[:,[0,2]]
+                rxClass = DC.Rx.Dipole_ky(rx[:,[0,2]], rx[:,[3,5]])
+            srcClass = DC.Src.Dipole([rxClass], M[0,:], N[-1,:])
+        SrcList.append(srcClass)
+    else:
+        print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
+
+
+    return SrcList
+
+
+def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
+    """
+        Write UBC GIF DCIP 2D or 3D observation file
+
+        Input:
+        :string     fileName -> including path where the file is written out
+        :DCsurvey   DC survey class object
+        :string     dtype ->  either '2D' | '3D'
+        :string     stype ->  either 'SURFACE' | 'GENERAL'
+
+        Output:
+        :param UBC2D-Data file
+        :return
+
+        Last edit: February 16th, 2016
+
+        @author: dominiquef
+
+    """
+    from SimPEG import mkvc
+
+    assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'"
+    assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
+
+    fid = open(fileName,'w')
+
+
+    if iptype!=0:
+        fid.write('IPTYPE=%i\n'%iptype)
+
+    else:
+        fid.write('! ' + stype + ' FORMAT\n')
+
+    count = 0
+
+    for ii in range(DCsurvey.nSrc):
+
+        tx = np.c_[DCsurvey.srcList[ii].loc]
+
+        rx = DCsurvey.srcList[ii].rxList[0].locs
+
+        nD = DCsurvey.srcList[ii].nD
+
+        M = rx[0]
+        N = rx[1]
+
+        # Adapt source-receiver location for dtype and stype
+        if dtype=='2D':
+
+            if stype == 'SIMPLE':
+
+                #fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
+                A = np.repeat(tx[0,0],M.shape[0],axis=0)
+                B = np.repeat(tx[0,1],M.shape[0],axis=0)
+                M = M[:,0]
+                N = N[:,0]
+
+                np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
+
+
+            else:
+
+                if stype == 'SURFACE':
+
+                    fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
+                    M = M[:,0]
+                    N = N[:,0]
+
+                if stype == 'GENERAL':
+
+                    # Flip sign for z-elevation to depth
+                    tx[2::2,:] = -tx[2::2,:]
+
+                    fid.writelines("%e " % ii for ii in mkvc(tx[::2,:]))
+                    M = M[:,0::2]
+                    N = N[:,0::2]
+
+                    # Flip sign for z-elevation to depth
+                    M[:,1::2] = -M[:,1::2]
+                    N[:,1::2] = -N[:,1::2]
+
+                fid.write('%i\n'% nD)
+                np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
+
+        if dtype=='3D':
+
+            if stype == 'SURFACE':
+
+                fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:]))
+                M = M[:,0:2]
+                N = N[:,0:2]
+
+            if stype == 'GENERAL':
+
+                fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
+
+            fid.write('%i\n'% nD)
+            np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
+            fid.write('\n')
+
+        count += nD
+
+    fid.close()

SimPEG/EM/Static/Utils/__init__.py

@@ -0,0 +1 @@
+from StaticUtils import *

SimPEG/EM/Static/__init__.py

@@ -0,0 +1,3 @@
+import DC
+import IP
+import SIP

SimPEG/EM/TDEM/SurveyTDEM.py

@@ -87,7 +87,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
     def getInitialFields(self, mesh):
         """Vertical magnetic dipole, magnetic vector potential"""
         if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"        
+            print ">> Step waveform: Non-zero initial condition"
             if mesh._meshType is 'CYL':
                 if mesh.isSymmetric:
                     MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
@@ -96,8 +96,8 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
             elif mesh._meshType is 'TENSOR':
                 MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
             else:
-                raise Exception('Unknown mesh for VMD') 
-            return {"b": mesh.edgeCurl*MVP}               
+                raise Exception('Unknown mesh for VMD')
+            return {"b": mesh.edgeCurl*MVP}
         elif self.waveformType == "GENERAL":
             print ">> General waveform: Zero initial condition"
             return {"b": np.zeros(mesh.nF)}
@@ -113,7 +113,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
         elif mesh._meshType is 'TENSOR':
             MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
         else:
-            raise Exception('Unknown mesh for VMD')            
+            raise Exception('Unknown mesh for VMD')
         return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
 
 
@@ -122,7 +122,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
         self.loc = loc
         self.radius = radius
         self.waveformType = waveformType
-        SrcTDEM.__init__(self,rxList)        
+        SrcTDEM.__init__(self,rxList)
 
     def getInitialFields(self, mesh):
         """Circular Loop, magnetic vector potential"""
@@ -153,7 +153,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
         elif mesh._meshType is 'TENSOR':
             MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
         else:
-            raise Exception('Unknown mesh for CircularLoop')            
+            raise Exception('Unknown mesh for CircularLoop')
         return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
 
 

SimPEG/EM/Utils/testingUtils.py

@@ -20,56 +20,61 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
     mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
 
     if useMu is True:
-        mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] 
+        mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))]
     else:
         mapping = Maps.ExpMap(mesh)
 
     x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid
     XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5))
-    Rx0 = EM.FDEM.Rx(XYZ, comp)
+    Rx0 = getattr(EM.FDEM.Rx, 'Point_' + comp[0])
+    if comp[2] == 'r':
+        real_or_imag = 'real'
+    elif comp[2] == 'i':
+        real_or_imag = 'imag'
+    rx0 = Rx0(XYZ, comp[1], 'imag')
 
     Src = []
 
     for SrcType in SrcList:
         if SrcType is 'MagDipole':
-            Src.append(EM.FDEM.Src.MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
+            Src.append(EM.FDEM.Src.MagDipole([rx0], freq=freq, loc=np.r_[0.,0.,0.]))
         elif SrcType is 'MagDipole_Bfield':
-            Src.append(EM.FDEM.Src.MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
+            Src.append(EM.FDEM.Src.MagDipole_Bfield([rx0], freq=freq, loc=np.r_[0.,0.,0.]))
         elif SrcType is 'CircularLoop':
-            Src.append(EM.FDEM.Src.CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
+            Src.append(EM.FDEM.Src.CircularLoop([rx0], freq=freq, loc=np.r_[0.,0.,0.]))
         elif SrcType is 'RawVec':
             if fdemType is 'e' or fdemType is 'b':
                 S_m = np.zeros(mesh.nF)
                 S_e = np.zeros(mesh.nE)
                 S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
                 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
-                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
+                Src.append(EM.FDEM.Src.RawVec([rx0], freq, S_m, mesh.getEdgeInnerProduct()*S_e))
 
             elif fdemType is 'h' or fdemType is 'j':
                 S_m = np.zeros(mesh.nE)
                 S_e = np.zeros(mesh.nF)
                 S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
                 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
-                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
+                Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*S_m, S_e))
 
     if verbose:
         print '  Fetching %s problem' % (fdemType)
 
     if fdemType == 'e':
         survey = EM.FDEM.Survey(Src)
-        prb = EM.FDEM.Problem_e(mesh, mapping=mapping)
+        prb = EM.FDEM.Problem3D_e(mesh, mapping=mapping)
 
     elif fdemType == 'b':
         survey = EM.FDEM.Survey(Src)
-        prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
+        prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping)
 
     elif fdemType == 'j':
         survey = EM.FDEM.Survey(Src)
-        prb = EM.FDEM.Problem_j(mesh, mapping=mapping)
+        prb = EM.FDEM.Problem3D_j(mesh, mapping=mapping)
 
     elif fdemType == 'h':
         survey = EM.FDEM.Survey(Src)
-        prb = EM.FDEM.Problem_h(mesh, mapping=mapping)
+        prb = EM.FDEM.Problem3D_h(mesh, mapping=mapping)
 
     else:
         raise NotImplementedError()
@@ -90,7 +95,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
     prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
     mesh = prb1.mesh
     print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
-    
+
     logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
     mu = np.ones(mesh.nC)*MU
 

SimPEG/EM/__init__.py

@@ -1,5 +1,6 @@
 import TDEM
 import FDEM
+import Static
 import Base
 import Analytics
 import Utils

SimPEG/Examples/EM_FDEM_1D_Inversion.py

@@ -42,8 +42,8 @@ def run(plotIt=True):
         ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
 
 
-    rxOffset=10. 
-    bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi')
+    rxOffset=10.
+    bzi = EM.FDEM.Rx.Point_b(np.array([[rxOffset, 0., 1e-3]]), orientation='z', component='imag')
 
     freqs = np.logspace(1,3,10)
     srcLoc = np.array([0., 0., 10.])
@@ -51,7 +51,7 @@ def run(plotIt=True):
     srcList = [EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z') for freq in freqs]
 
     survey = EM.FDEM.Survey(srcList)
-    prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
+    prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping)
 
     try:
         from pymatsolver import MumpsSolver

SimPEG/Examples/EM_Schenkel_Morrison_Casing.py

@@ -215,7 +215,7 @@ def run(plotIt=True):
     # ------------ Problem and Survey ---------------
     survey = FDEM.Survey(sg_p + dg_p)
     mapping = [('sigma', Maps.IdentityMap(mesh))]
-    problem = FDEM.Problem_h(mesh, mapping=mapping)
+    problem = FDEM.Problem3D_h(mesh, mapping=mapping)
     problem.pair(survey)
 
     # ------------- Solve ---------------------------

SimPEG/MT/BaseMT.py

@@ -1,5 +1,5 @@
 from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
-from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem
+from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
 from SurveyMT import Survey, Data
 from FieldsMT import BaseMTFields
 

SimPEG/Mesh/CylMesh.py

@@ -330,7 +330,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
                 raise NotImplementedError('wrapping in the averaging is not yet implemented')
         return self._aveF2CCV
 
-    def getInterpolationMatCartMesh(self, Mrect, locType='CC'):
+    def getInterpolationMatCartMesh(self, Mrect, locType='CC', locTypeTo=None):
         """
             Takes a cartesian mesh and returns a projection to translate onto the cartesian grid.
         """
@@ -338,19 +338,22 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes"
 
 
+        if locTypeTo is None:
+            locTypeTo = locType
+
         if locType == 'F':
             # do this three times for each component
-            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx')
-            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy')
-            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz')
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx', locTypeTo=locTypeTo+'x')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy', locTypeTo=locTypeTo+'y')
+            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz', locTypeTo=locTypeTo+'z')
             return sp.vstack((X,Y,Z))
         if locType == 'E':
-            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex')
-            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey')
-            Z = spzeros(Mrect.nEz, self.nE)
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex', locTypeTo=locTypeTo+'x')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey', locTypeTo=locTypeTo+'y')
+            Z = spzeros(getattr(Mrect, 'n' + locTypeTo + 'z'), self.nE)
             return sp.vstack((X,Y,Z))
 
-        grid = getattr(Mrect, 'grid' + locType)
+        grid = getattr(Mrect, 'grid' + locTypeTo)
         # This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice
         theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2
         theta[theta < 0] += np.pi*2.0
@@ -366,7 +369,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
                         'Ex': Mrect.tangents[:Mrect.nEx,:],
                         'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:],
                         'Ez': Mrect.tangents[-Mrect.nEz:,:],
-                    }[locType]
+                    }[locTypeTo]
             if 'F' in locType:
                 normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)]
                 proj = ( normals * dotMe ).sum(axis=1)

SimPEG/Mesh/DiffOperators.py

@@ -584,7 +584,67 @@ class DiffOperators(object):
 
         return Pbc, Pin, Pout
 
+    def getBCProjWF_simple(self, discretization='CC'):
+        """
 
+        The weak form boundary condition projection matrices
+        when mixed boundary condition is used
+
+
+        """
+
+        if discretization is not 'CC':
+            raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
+
+        def projBC(n):
+            ij = ([0,n], [0,1])
+            vals = [0,0]
+            vals[0] = 1
+            vals[1] = 1
+            return sp.csr_matrix((vals, ij), shape=(n+1,2))
+
+        def projDirichlet(n, bc):
+            bc = checkBC(bc)
+            ij = ([0,n], [0,1])
+            vals = [0,0]
+            if(bc[0] == 'dirichlet'):
+                vals[0] = -1
+            if(bc[1] == 'dirichlet'):
+                vals[1] = 1
+            return sp.csr_matrix((vals, ij), shape=(n+1,2))
+
+        BC = [['dirichlet','dirichlet'],['dirichlet','dirichlet'],['dirichlet','dirichlet']]
+        n = self.vnC
+        indF = self.faceBoundaryInd
+        if(self.dim == 1):
+            Pbc = projDirichlet(n[0], BC[0])
+            B = projBC(n[0])
+            indF = indF[0] | indF[1]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+        elif(self.dim == 2):
+            Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
+            Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
+            Pbc = sp.block_diag((Pbc1, Pbc2), format="csr")
+            B1 = sp.kron(speye(n[1]), projBC(n[0]))
+            B2 = sp.kron(projBC(n[1]), speye(n[0]))
+            B = sp.block_diag((B1, B2), format="csr")
+            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+        elif(self.dim == 3):
+            Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
+            Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
+            Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
+            Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
+            B1 = kron3(speye(n[2]), speye(n[1]), projBC(n[0]))
+            B2 = kron3(speye(n[2]), projBC(n[1]), speye(n[0]))
+            B3 = kron3(projBC(n[2]), speye(n[1]), speye(n[0]))
+            B = sp.block_diag((B1, B2, B3), format="csr")
+            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+        return Pbc, B.T
     # --------------- Averaging ---------------------
 
     @property

SimPEG/Mesh/MeshIO.py

@@ -205,19 +205,30 @@ class TensorMeshIO(object):
             :param simpeg.Mesh.TensorMesh mesh: The mesh
 
         """
-        assert mesh.dim == 3
-        s = ''
-        s += '%i %i %i\n' %tuple(mesh.vnC)
-        origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated.
-        origin.dtype = float
+        if mesh.dim ==3:
+            s = ''
+            s += '%i %i %i\n' %tuple(mesh.vnC)
+            origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated.
+            origin.dtype = float
 
-        s += '%.2f %.2f %.2f\n' %tuple(origin)
-        s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
-        s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
-        s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1])
-        f = open(fileName, 'w')
-        f.write(s)
-        f.close()
+            s += '%.2f %.2f %.2f\n' %tuple(origin)
+            s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
+            s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
+            s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1])
+            f = open(fileName, 'w')
+            f.write(s)
+            f.close()
+
+        elif mesh.dim==2:
+            fid = open(fileName,'w')
+            fid.write('%i\n'% mesh.nCx)
+            fid.write('%f %f 1\n'% (mesh.vectorNx[0],mesh.vectorNx[1]))
+            np.savetxt(fid, np.c_[mesh.vectorNx[2:],np.ones(mesh.nCx-1)], fmt='\t %e %i',delimiter=' ',newline='\n')
+            fid.write('\n')
+            fid.write('%i\n'% mesh.nCy)
+            fid.write('%f %f 1\n'%( 0,mesh.hy[-1]))
+            np.savetxt(fid, np.c_[np.cumsum(mesh.hy[-2::-1])+mesh.hy[-1],np.ones(mesh.nCy-1)], fmt='\t %e %i',delimiter=' ',newline='\n')
+            fid.close()
 
         if models is None: return
         assert type(models) is dict, 'models must be a dict'

SimPEG/Mesh/View.py

@@ -218,7 +218,7 @@ class TensorView(object):
             return out
         viewOpts = ['real','imag','abs','vec']
         normalOpts = ['X', 'Y', 'Z']
-        vTypeOpts = ['CC', 'CCv','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez']
+        vTypeOpts = ['CC', 'CCv','N','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez']
 
         # Some user error checking
         assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)

SimPEG/Regularization.py

@@ -311,6 +311,9 @@ class BaseRegularization(object):
             tmp = indActive
             indActive = np.zeros(mesh.nC, dtype=bool)
             indActive[tmp] = True
+        if indActive is not None and mapping is None:
+            mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
+
         self.regmesh = RegularizationMesh(mesh,indActive)
         self.mapping = mapping or self.mapPair(mesh)
         self.mapping._assertMatchesPair(self.mapPair)

docs/examples/DC_Forward_PseudoSection.rst

@@ -22,7 +22,7 @@ radi    = Radius of spheres [r1,r2]
 param   = Conductivity of background and two spheres [m0,m1,m2]
 stype   = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
 dtype   = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
-Created by @fourndo on Mon Feb 01 19:28:06 2016
+Created by @fourndo
 
 
 

setup.py

@@ -5,16 +5,17 @@ SimPEG is a python package for simulation and gradient based
 parameter estimation in the context of geophysical applications.
 """
 
-import numpy as np
-
 import os
 import sys
 import subprocess
 
 from distutils.core import setup
+from distutils.command.build_ext import build_ext
 from setuptools import find_packages
 from distutils.extension import Extension
 
+
+
 CLASSIFIERS = [
 'Development Status :: 4 - Beta',
 'Intended Audience :: Developers',
@@ -51,11 +52,16 @@ if args.count("build_ext") > 0 and args.count("--inplace") == 0:
 try:
     from Cython.Build import cythonize
     from Cython.Distutils import build_ext
-    cythonKwargs = dict(cmdclass={'build_ext': build_ext})
     USE_CYTHON = True
 except Exception, e:
     USE_CYTHON = False
-    cythonKwargs = dict()
+
+class NumpyBuild(build_ext):
+    def finalize_options(self):
+        build_ext.finalize_options(self)
+        __builtins__.__NUMPY_SETUP__ = False
+        import numpy
+        self.include_dirs.append(numpy.get_include())
 
 ext = '.pyx' if USE_CYTHON else '.c'
 
@@ -94,8 +100,8 @@ setup(
     classifiers=CLASSIFIERS,
     platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"],
     use_2to3 = False,
-    include_dirs=[np.get_include()],
+    cmdclass={'build_ext':NumpyBuild},
+    setup_requires=['numpy'],
     ext_modules = extensions,
     scripts=scripts,
-    **cythonKwargs
 )

tests/base/test_regularization.py

@@ -65,10 +65,8 @@ class RegularizationTests(unittest.TestCase):
                     elif mesh.dim == 3:
                         indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
 
-                    mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
-
                     for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers
-                        reg = r(mesh, mapping=mapping, indActive=indAct)
+                        reg = r(mesh, indActive=indAct)
                         m = np.random.rand(mesh.nC)[indAct]
                         reg.mref = np.ones_like(m)*np.mean(m)
 

tests/em/fdem/forward/test_FDEM_analytics.py

@@ -28,12 +28,12 @@ class FDEM_analyticTests(unittest.TestCase):
 
         x = np.linspace(-10,10,5)
         XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
-        rxList = EM.FDEM.Rx(XYZ, 'exi')
+        rxList = EM.FDEM.Rx.Point_e(XYZ, orientation='x', component='imag')
         Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
 
         survey = EM.FDEM.Survey([Src0])
 
-        prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
+        prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping)
         prb.pair(survey)
 
         try:
@@ -125,8 +125,8 @@ class FDEM_analyticTests(unittest.TestCase):
 
         mapping = [('sigma', Maps.IdentityMap(mesh)),('mu', Maps.IdentityMap(mesh))]
 
-        prbe = EM.FDEM.Problem_h(mesh, mapping=mapping)
-        prbm = EM.FDEM.Problem_e(mesh, mapping=mapping)
+        prbe = EM.FDEM.Problem3D_h(mesh, mapping=mapping)
+        prbm = EM.FDEM.Problem3D_e(mesh, mapping=mapping)
 
         prbe.pair(surveye) # pair problem and survey
         prbm.pair(surveym)

tests/em/fdem/forward/test_FDEM_forwardHB.py

@@ -12,7 +12,7 @@ testBH = True
 verbose = False
 
 TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.)
-#TODO: choose better testing parameters to lower this 
+#TODO: choose better testing parameters to lower this
 
 SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop']
 
@@ -125,4 +125,4 @@ class FDEM_CrossCheck(unittest.TestCase):
             self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB))
 
 if __name__ == '__main__':
-    unittest.main()
+    unittest.main()

tests/em/static/__init__.py

@@ -0,0 +1,12 @@
+import os
+import glob
+import unittest
+
+if __name__ == '__main__':
+    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)

tests/em/static/test_DC_2D_analytic.py

@@ -0,0 +1,69 @@
+import unittest
+from SimPEG import Mesh, Utils, EM, Maps, np
+import SimPEG.EM.Static.DC as DC
+
+class DCProblemAnalyticTests(unittest.TestCase):
+
+    def setUp(self):
+
+        cs = 12.5
+        hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)]
+        hy = [(cs,7, -1.3),(cs,20)]
+        mesh = Mesh.TensorMesh([hx, hy],x0="CN")
+        sighalf = 1e-2
+        sigma = np.ones(mesh.nC)*sighalf
+        x = np.linspace(-135, 250., 20)
+        M = Utils.ndgrid(x-12.5, np.r_[0.])
+        N = Utils.ndgrid(x+12.5, np.r_[0.])
+        A0loc = np.r_[-150, 0.]
+        A1loc = np.r_[-130, 0.]
+        rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
+        data_anal = EM.Analytics.DCAnalyticHalf(np.r_[A0loc, 0.], rxloc, sighalf, earth_type="halfspace")
+
+        rx = DC.Rx.Dipole_ky(M, N)
+        src0 = DC.Src.Pole([rx], A0loc)
+        survey = DC.Survey_ky([src0])
+
+        self.survey = survey
+        self.mesh = mesh
+        self.sigma = sigma
+        self.data_anal = data_anal
+
+        try:
+            from pymatsolver import MumpsSolver
+            self.Solver = MumpsSolver
+        except ImportError, e:
+            self.Solver = SolverLU
+
+    def test_Problem3D_N(self):
+
+        problem = DC.Problem2D_N(self.mesh)
+        problem.Solver = self.Solver
+        problem.pair(self.survey)
+        data = self.survey.dpred(self.sigma)
+        err= np.linalg.norm((data-self.data_anal)/self.data_anal)**2 / self.data_anal.size
+        if err < 0.05:
+            passed = True
+            print ">> DC analytic test for Problem3D_N is passed"
+        else:
+            passed = False
+            print ">> DC analytic test for Problem3D_N is failed"
+        self.assertTrue(passed)
+
+    def test_Problem3D_CC(self):
+        problem = DC.Problem2D_CC(self.mesh)
+        problem.Solver = self.Solver
+        problem.pair(self.survey)
+        data = self.survey.dpred(self.sigma)
+        err= np.linalg.norm((data-self.data_anal)/self.data_anal)**2 / self.data_anal.size
+        if err < 0.05:
+            passed = True
+            print ">> DC analytic test for Problem3D_CC is passed"
+        else:
+            passed = False
+            print ">> DC analytic test for Problem3D_CC is failed"
+        self.assertTrue(passed)
+
+if __name__ == '__main__':
+    unittest.main()
+

tests/em/static/test_DC_2D_jvecjtvecadj.py

@@ -0,0 +1,127 @@
+import unittest
+from SimPEG import *
+import SimPEG.EM.Static.DC as DC
+
+
+class DCProblem_2DTestsCC(unittest.TestCase):
+
+    def setUp(self):
+
+        cs = 12.5
+        hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)]
+        hy = [(cs,7, -1.3),(cs,20)]
+        mesh = Mesh.TensorMesh([hx, hy],x0="CN")
+        x = np.linspace(-135, 250., 20)
+        M = Utils.ndgrid(x-12.5, np.r_[0.])
+        N = Utils.ndgrid(x+12.5, np.r_[0.])
+        A0loc = np.r_[-150, 0.]
+        A1loc = np.r_[-130, 0.]
+        rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
+        rx = DC.Rx.Dipole_ky(M, N)
+        src0 = DC.Src.Pole([rx], A0loc)
+        src1 = DC.Src.Pole([rx], A1loc)
+        survey = DC.Survey_ky([src0, src1])
+        problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
+        problem.pair(survey)
+
+        mSynth = np.ones(mesh.nC)*1.
+        survey.makeSyntheticData(mSynth)
+
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv    = inv
+        self.reg    = reg
+        self.p      = problem
+        self.mesh   = mesh
+        self.m0     = mSynth
+        self.survey = survey
+        self.dmis   = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-10
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+class DCProblemTestsN(unittest.TestCase):
+
+    def setUp(self):
+
+        cs = 12.5
+        hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)]
+        hy = [(cs,7, -1.3),(cs,20)]
+        mesh = Mesh.TensorMesh([hx, hy],x0="CN")
+        x = np.linspace(-135, 250., 20)
+        M = Utils.ndgrid(x-12.5, np.r_[0.])
+        N = Utils.ndgrid(x+12.5, np.r_[0.])
+        A0loc = np.r_[-150, 0.]
+        A1loc = np.r_[-130, 0.]
+        rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
+        rx = DC.Rx.Dipole_ky(M, N)
+        src0 = DC.Src.Pole([rx], A0loc)
+        src1 = DC.Src.Pole([rx], A1loc)
+        survey = DC.Survey_ky([src0, src1])
+        problem = DC.Problem2D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
+        problem.pair(survey)
+
+        mSynth = np.ones(mesh.nC)*1.
+        survey.makeSyntheticData(mSynth)
+
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p =     problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-8
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+if __name__ == '__main__':
+    unittest.main()

tests/em/static/test_DC_analytic.py

@@ -0,0 +1,71 @@
+import unittest
+from SimPEG import Mesh, Utils, EM, Maps, np
+import SimPEG.EM.Static.DC as DC
+
+class DCProblemAnalyticTests(unittest.TestCase):
+
+    def setUp(self):
+
+        cs = 25.
+        hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
+        hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
+        hz = [(cs,7, -1.3),(cs,20)]
+        mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
+        sigma = np.ones(mesh.nC)*1e-2
+
+        x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
+        y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
+        Aloc = np.r_[-200., 0., 0.]
+        Bloc = np.r_[200., 0., 0.]
+        M = Utils.ndgrid(x-25.,y, np.r_[0.])
+        N = Utils.ndgrid(x+25.,y, np.r_[0.])
+        phiA = EM.Analytics.DCAnalyticHalf(Aloc, [M,N], 1e-2, earth_type="halfspace")
+        phiB = EM.Analytics.DCAnalyticHalf(Bloc, [M,N], 1e-2, earth_type="halfspace")
+        data_anal = phiA-phiB
+
+        rx = DC.Rx.Dipole(M, N)
+        src = DC.Src.Dipole([rx], Aloc, Bloc)
+        survey = DC.Survey([src])
+
+        self.survey = survey
+        self.mesh = mesh
+        self.sigma = sigma
+        self.data_anal = data_anal
+
+        try:
+            from pymatsolver import MumpsSolver
+            self.Solver = MumpsSolver
+        except ImportError, e:
+            self.Solver = SolverLU
+
+    def test_Problem3D_N(self):
+        problem = DC.Problem3D_N(self.mesh)
+        problem.Solver = self.Solver
+        problem.pair(self.survey)
+        data = self.survey.dpred(self.sigma)
+        err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal)
+        if err < 0.2:
+            passed = True
+            print ">> DC analytic test for Problem3D_N is passed"
+        else:
+            passed = False
+            print ">> DC analytic test for Problem3D_N is failed"
+        self.assertTrue(passed)
+
+    def test_Problem3D_CC(self):
+        problem = DC.Problem3D_CC(self.mesh)
+        problem.Solver = self.Solver
+        problem.pair(self.survey)
+        data = self.survey.dpred(self.sigma)
+        err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal)
+        if err < 0.2:
+            passed = True
+            print ">> DC analytic test for Problem3D_CC is passed"
+        else:
+            passed = False
+            print ">> DC analytic test for Problem3D_CC is failed"
+        self.assertTrue(passed)
+
+if __name__ == '__main__':
+    unittest.main()
+

tests/em/static/test_DC_jvecjtvecadj.py

@@ -0,0 +1,127 @@
+import unittest
+from SimPEG import *
+import SimPEG.EM.Static.DC as DC
+
+
+class DCProblemTestsCC(unittest.TestCase):
+
+    def setUp(self):
+
+        aSpacing=2.5
+        nElecs=5
+
+        surveySize = nElecs*aSpacing - aSpacing
+        cs = surveySize/nElecs/4
+
+        mesh = Mesh.TensorMesh([
+                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
+                [(cs,3, -1.3),(cs,3,1.3)],
+                # [(cs,5, -1.3),(cs,10)]
+            ],'CN')
+
+        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
+        survey = DC.Survey(srcList)
+        problem = DC.Problem3D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
+        problem.pair(survey)
+
+        mSynth = np.ones(mesh.nC)
+        survey.makeSyntheticData(mSynth)
+
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p =     problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-10
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+class DCProblemTestsN(unittest.TestCase):
+
+    def setUp(self):
+
+        aSpacing=2.5
+        nElecs=10
+
+        surveySize = nElecs*aSpacing - aSpacing
+        cs = surveySize/nElecs/4
+
+        mesh = Mesh.TensorMesh([
+                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
+                [(cs,3, -1.3),(cs,3,1.3)],
+        #         [(cs,5, -1.3),(cs,10)]
+            ],'CN')
+
+        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
+        survey = DC.Survey(srcList)
+        problem = DC.Problem3D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
+        problem.pair(survey)
+
+        mSynth = np.ones(mesh.nC)
+        survey.makeSyntheticData(mSynth)
+
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p =     problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-8
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+if __name__ == '__main__':
+    unittest.main()

tests/em/static/test_IP_fwd.py

@@ -0,0 +1,96 @@
+import unittest
+from SimPEG import Mesh, Utils, EM, Maps, np
+import SimPEG.EM.Static.DC as DC
+import SimPEG.EM.Static.IP as IP
+
+class IPProblemAnalyticTests(unittest.TestCase):
+
+    def setUp(self):
+
+        cs = 12.5
+        npad=2
+        hx = [(cs,npad, -1.3),(cs,21),(cs,npad, 1.3)]
+        hy = [(cs,npad, -1.3),(cs,21),(cs,npad, 1.3)]
+        hz = [(cs,npad, -1.3),(cs,20)]
+        mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
+
+        x = mesh.vectorCCx[(mesh.vectorCCx>-80.)&(mesh.vectorCCx<80.)]
+        y = mesh.vectorCCx[(mesh.vectorCCy>-80.)&(mesh.vectorCCy<80.)]
+        Aloc = np.r_[-100., 0., 0.]
+        Bloc = np.r_[100., 0., 0.]
+        M = Utils.ndgrid(x-12.5,y, np.r_[0.])
+        N = Utils.ndgrid(x+12.5,y, np.r_[0.])
+        radius = 50.
+        xc = np.r_[0., 0., -100]
+        blkind = Utils.ModelBuilder.getIndicesSphere(xc, radius, mesh.gridCC)
+        sigmaInf = np.ones(mesh.nC)*1e-2
+        eta = np.zeros(mesh.nC)
+        eta[blkind] = 0.1
+        sigma0 = sigmaInf*(1.-eta)
+
+        rx = DC.Rx.Dipole(M, N)
+        src = DC.Src.Dipole([rx], Aloc, Bloc)
+        surveyDC = DC.Survey([src])
+
+        self.surveyDC = surveyDC
+        self.mesh = mesh
+        self.sigmaInf = sigmaInf
+        self.sigma0 = sigma0
+        self.src = src
+        self.eta = eta
+
+        try:
+            from pymatsolver import MumpsSolver
+            self.Solver = MumpsSolver
+        except ImportError, e:
+            self.Solver = SolverLU
+
+    def test_Problem3D_N(self):
+
+        problemDC = DC.Problem3D_N(self.mesh)
+        problemDC.Solver = self.Solver
+        problemDC.pair(self.surveyDC)
+        data0 = self.surveyDC.dpred(self.sigma0)
+        finf = problemDC.fields(self.sigmaInf)
+        datainf = self.surveyDC.dpred(self.sigmaInf, f=finf)
+        problemIP = IP.Problem3D_N(self.mesh, sigma=self.sigmaInf, Ainv=problemDC.Ainv, f=finf)
+        problemIP.Solver = self.Solver
+        surveyIP = IP.Survey([self.src])
+        problemIP.pair(surveyIP)
+        data_full = data0 - datainf
+        data = surveyIP.dpred(self.eta)
+        err= np.linalg.norm((data-data_full)/data_full)**2 / data_full.size
+        if err < 0.05:
+            passed = True
+            print ">> IP forward test for Problem3D_N is passed"
+        else:
+            passed = False
+            print ">> IP forward test for Problem3D_N is failed"
+        self.assertTrue(passed)
+
+    def test_Problem3D_CC(self):
+
+        problemDC = DC.Problem3D_CC(self.mesh)
+        problemDC.Solver = self.Solver
+        problemDC.pair(self.surveyDC)
+        data0 = self.surveyDC.dpred(self.sigma0)
+        finf = problemDC.fields(self.sigmaInf)
+        datainf = self.surveyDC.dpred(self.sigmaInf, f=finf)
+        problemIP = IP.Problem3D_CC(self.mesh, rho=1./self.sigmaInf, Ainv=problemDC.Ainv, f=finf)
+        problemIP.Solver = self.Solver
+        surveyIP = IP.Survey([self.src])
+        problemIP.pair(surveyIP)
+        data_full = data0 - datainf
+        data = surveyIP.dpred(self.eta)
+        err= np.linalg.norm((data-data_full)/data_full)**2 / data_full.size
+        if err < 0.05:
+            passed = True
+            print ">> IP forward test for Problem3D_CC is passed"
+        else:
+            passed = False
+            print ">> IP forward test for Problem3D_CC is failed"
+        self.assertTrue(passed)
+
+if __name__ == '__main__':
+    unittest.main()
+

tests/em/static/test_IP_jvecjtvecadj.py

@@ -0,0 +1,126 @@
+import unittest
+from SimPEG import *
+import SimPEG.EM.Static.DC as DC
+import SimPEG.EM.Static.IP as IP
+
+
+class IPProblemTestsCC(unittest.TestCase):
+
+    def setUp(self):
+
+        aSpacing=2.5
+        nElecs=5
+
+        surveySize = nElecs*aSpacing - aSpacing
+        cs = surveySize/nElecs/4
+
+        mesh = Mesh.TensorMesh([
+                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
+                [(cs,3, -1.3),(cs,3,1.3)],
+                # [(cs,5, -1.3),(cs,10)]
+            ],'CN')
+
+        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
+        survey = IP.Survey(srcList)
+        sigma = np.ones(mesh.nC)
+        problem = IP.Problem3D_CC(mesh, rho=1./sigma)
+        problem.pair(survey)
+        mSynth = np.ones(mesh.nC)*0.1
+        survey.makeSyntheticData(mSynth)
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p =     problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-10
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+class IPProblemTestsN(unittest.TestCase):
+
+    def setUp(self):
+
+        aSpacing=2.5
+        nElecs=5
+
+        surveySize = nElecs*aSpacing - aSpacing
+        cs = surveySize/nElecs/4
+
+        mesh = Mesh.TensorMesh([
+                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
+                [(cs,3, -1.3),(cs,3,1.3)],
+                # [(cs,5, -1.3),(cs,10)]
+            ],'CN')
+
+        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
+        survey = IP.Survey(srcList)
+        sigma = np.ones(mesh.nC)
+        problem = IP.Problem3D_N(mesh, sigma=sigma)
+        problem.pair(survey)
+        mSynth = np.ones(mesh.nC)*0.1
+        survey.makeSyntheticData(mSynth)
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p =     problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-8
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+if __name__ == '__main__':
+    unittest.main()

tests/em/static/test_SIP_jvecjtvecadj.py

@@ -0,0 +1,232 @@
+import unittest
+from SimPEG import *
+import SimPEG
+from SimPEG import Mesh, Utils, EM, Maps, np, Survey
+from SimPEG.EM.Static import SIP, DC, IP
+from pymatsolver import MumpsSolver
+
+
+class IPProblemTestsCC(unittest.TestCase):
+
+    def setUp(self):
+
+        cs = 25.
+        hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
+        hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
+        hz = [(cs,0, -1.3),(cs,20)]
+        mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
+        blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
+        blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
+        sigma = np.ones(mesh.nC)*1e-2
+        eta = np.zeros(mesh.nC)
+        tau = np.ones_like(sigma)*1.
+        eta[blkind0] = 0.1
+        eta[blkind1] = 0.1
+        tau[blkind0] = 0.1
+        tau[blkind1] = 0.01
+
+        x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
+        y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
+        Aloc = np.r_[-200., 0., 0.]
+        Bloc = np.r_[200., 0., 0.]
+        M = Utils.ndgrid(x-25.,y, np.r_[0.])
+        N = Utils.ndgrid(x+25.,y, np.r_[0.])
+
+        times = np.arange(10)*1e-3 + 1e-3
+        rx = SIP.Rx.Dipole(M, N, times)
+        src = SIP.Src.Dipole([rx], Aloc, Bloc)
+        survey = SIP.Survey([src])
+        colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))]
+        problem = SIP.Problem3D_CC(mesh, rho=1./sigma, mapping=colemap)
+        problem.Solver = MumpsSolver
+        problem.pair(survey)
+        mSynth = np.r_[eta, 1./tau]
+        survey.makeSyntheticData(mSynth)
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p =     problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC*2)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-10
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+class IPProblemTestsN(unittest.TestCase):
+
+    def setUp(self):
+
+        cs = 25.
+        hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
+        hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
+        hz = [(cs,0, -1.3),(cs,20)]
+        mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
+        blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
+        blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
+        sigma = np.ones(mesh.nC)*1e-2
+        eta = np.zeros(mesh.nC)
+        tau = np.ones_like(sigma)*1.
+        eta[blkind0] = 0.1
+        eta[blkind1] = 0.1
+        tau[blkind0] = 0.1
+        tau[blkind1] = 0.01
+
+        x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
+        y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
+        Aloc = np.r_[-200., 0., 0.]
+        Bloc = np.r_[200., 0., 0.]
+        M = Utils.ndgrid(x-25.,y, np.r_[0.])
+        N = Utils.ndgrid(x+25.,y, np.r_[0.])
+
+        times = np.arange(10)*1e-3 + 1e-3
+        rx = SIP.Rx.Dipole(M, N, times)
+        src = SIP.Src.Dipole([rx], Aloc, Bloc)
+        survey = SIP.Survey([src])
+        colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))]
+        problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap)
+        problem.Solver = MumpsSolver
+        problem.pair(survey)
+        mSynth = np.r_[eta, 1./tau]
+        survey.makeSyntheticData(mSynth)
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p =     problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC*2)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-8
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+class IPProblemTestsN_air(unittest.TestCase):
+
+    def setUp(self):
+
+        cs = 25.
+        hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
+        hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
+        hz = [(cs,0, -1.3),(cs,20),(cs,0, 1.3)]
+        mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCC")
+        blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
+        blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
+        sigma = np.ones(mesh.nC)*1e-2
+        airind = mesh.gridCC[:,2]>0.
+        sigma[airind] = 1e-8
+        eta = np.zeros(mesh.nC)
+        tau = np.ones_like(sigma)*1.
+        eta[blkind0] = 0.1
+        eta[blkind1] = 0.1
+        tau[blkind0] = 0.1
+        tau[blkind1] = 0.01
+
+        actmapeta = Maps.InjectActiveCells(mesh, ~airind, 0.)
+        actmaptau = Maps.InjectActiveCells(mesh, ~airind, 1.)
+
+        x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
+        y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
+        Aloc = np.r_[-200., 0., 0.]
+        Bloc = np.r_[200., 0., 0.]
+        M = Utils.ndgrid(x-25.,y, np.r_[0.])
+        N = Utils.ndgrid(x+25.,y, np.r_[0.])
+
+        times = np.arange(10)*1e-3 + 1e-3
+        rx = SIP.Rx.Dipole(M, N, times)
+        src = SIP.Src.Dipole([rx], Aloc, Bloc)
+        survey = SIP.Survey([src])
+        colemap = [("eta", Maps.IdentityMap(mesh)*actmapeta), ("taui", Maps.IdentityMap(mesh)*actmaptau)]
+        problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap)
+        problem.Solver = MumpsSolver
+        problem.pair(survey)
+        mSynth = np.r_[eta[~airind], 1./tau[~airind]]
+        survey.makeSyntheticData(mSynth)
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        regmap = Maps.IdentityMap(nP=int(mSynth[~airind].size*2))
+        reg = SIP.MultiRegularization(mesh, mapping=regmap, nModels=2, indActive=~airind)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p =     problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-8
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
+        self.assertTrue(passed)
+
+if __name__ == '__main__':
+    unittest.main()

tests/mesh/test_Mixed_boundaryPoisson.py

@@ -0,0 +1,411 @@
+import numpy as np
+import scipy.sparse as sp
+import unittest
+import matplotlib.pyplot as plt
+from SimPEG import *
+
+MESHTYPES = ['uniformTensorMesh']
+
+def getxBCyBC_CC(mesh, alpha, beta, gamma):
+# def getxBCyBC(mesh, alpha, beta, gamma):
+    """
+    This is a subfunction generating mixed-boundary condition:
+
+    .. math::
+
+        \nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q
+
+        \rho \vec{j} = -\nabla \phi \phi
+
+        \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega
+
+        xBC = f_1(\alpha, \beta, \gamma)
+        yBC = f(\alpha, \beta, \gamma)
+
+    Computes xBC and yBC for cell-centered discretizations
+    """
+    if mesh.dim == 1: #1D
+        if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2):
+            raise Exception("Lenght of list, alpha should be 2")
+        fCCxm,fCCxp = mesh.cellBoundaryInd
+        nBC = fCCxm.sum()+fCCxp.sum()
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+
+        # h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+        h_xm, h_xp = mesh.hx[0], mesh.hx[-1]
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+
+        xBC = np.r_[xBC_xm, xBC_xp]
+        yBC = np.r_[yBC_xm, yBC_xp]
+
+    elif mesh.dim == 2: #2D
+        if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4):
+            raise Exception("Lenght of list, alpha should be 4")
+
+        fxm,fxp,fym,fyp = mesh.faceBoundaryInd
+        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
+        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
+
+        # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
+        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
+
+        h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
+        h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
+        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
+        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
+        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+        xBC_ym = 0.5*a_ym
+        xBC_yp = 0.5*a_yp/b_yp
+        yBC_ym = 0.5*(1.-b_ym)
+        yBC_yp = 0.5*(1.-1./b_yp)
+
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
+
+        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
+        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
+        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
+        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
+
+        xBC = np.r_[xBC_x, xBC_y]
+        yBC = np.r_[yBC_x, yBC_y]
+
+    elif mesh.dim == 3: #3D
+        if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6):
+            raise Exception("Lenght of list, alpha should be 6")
+        # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd
+        fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd
+        nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
+        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
+        alpha_zm, beta_zm, gamma_zm = alpha[4], beta[4], gamma[4]
+        alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5]
+
+        # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
+        # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
+        # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2]
+
+        h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
+        h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
+        h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp)
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
+        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
+        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
+        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
+
+        a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm)
+        b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm)
+        a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp)
+        b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+        xBC_ym = 0.5*a_ym
+        xBC_yp = 0.5*a_yp/b_yp
+        yBC_ym = 0.5*(1.-b_ym)
+        yBC_yp = 0.5*(1.-1./b_yp)
+        xBC_zm = 0.5*a_zm
+        xBC_zp = 0.5*a_zp/b_zp
+        yBC_zm = 0.5*(1.-b_zm)
+        yBC_zp = 0.5*(1.-1./b_zp)
+
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
+        sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]])
+
+        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
+        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
+        xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz]
+
+        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
+        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
+        yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz]
+
+        xBC = np.r_[xBC_x, xBC_y, xBC_z]
+        yBC = np.r_[yBC_x, yBC_y, yBC_z]
+
+    return xBC, yBC
+
+class Test1D_InhomogeneousMixed(Tests.OrderTest):
+    name = "1D - Mixed"
+    meshTypes = MESHTYPES
+    meshDimension = 1
+    expectedOrders = 2
+    meshSizes = [4, 8, 16, 32]
+
+    def getError(self):
+        #Test function
+        phi_fun = lambda x: np.cos(np.pi*x)
+        j_fun = lambda x: np.pi*np.sin(np.pi*x)
+        phi_deriv = lambda x: -j_fun(x)
+        q_fun = lambda x: (np.pi**2)*np.cos(np.pi*x)
+
+        xc_ana = phi_fun(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        j_ana = j_fun(self.M.gridFx)
+
+        # Get boundary locations
+        vecN = self.M.vectorNx
+        vecC = self.M.vectorCCx
+
+        # Setup Mixed B.C (alpha, beta, gamma)
+        alpha_xm, alpha_xp = 1., 1.
+        beta_xm, beta_xp =  1., 1.
+        alpha = np.r_[alpha_xm, alpha_xp]
+        beta = np.r_[beta_xm, beta_xp]
+        vecN = self.M.vectorNx
+        vecC = self.M.vectorCCx
+        phi_bc = phi_fun(vecN[[0,-1]])
+        phi_deriv_bc = phi_deriv(vecN[[0,-1]])
+        gamma = alpha*phi_bc + beta*phi_deriv_bc
+        x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
+
+
+        sigma = np.ones(self.M.nC)
+        Mfrho = self.M.getFaceInnerProduct(1./sigma)
+        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
+        V = Utils.sdiag(self.M.vol)
+        Div = V*self.M.faceDiv
+        P_BC, B = self.M.getBCProjWF_simple()
+        q = q_fun(self.M.gridCC)
+        M = B*self.M.aveCC2F
+        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
+        # Mrhoj = D.T V phi + P_BC*Utils.sdiag(y_BC)*M phi - P_BC*x_BC
+        rhs = V*q + Div*MfrhoI*P_BC*x_BC
+        A = Div*MfrhoI*G
+
+        if self.myTest == 'xc':
+            #TODO: fix the null space
+            Ainv = Solver(A)
+            xc = Ainv*rhs
+            err = np.linalg.norm((xc-xc_ana), np.inf)
+        else:
+            NotImplementedError
+        return err
+
+
+    def test_order(self):
+        print "==== Testing Mixed boudary conduction for CC-problem ===="
+        self.name = "1D"
+        self.myTest = 'xc'
+        self.orderTest()
+
+class Test2D_InhomogeneousMixed(Tests.OrderTest):
+    name = "2D - Mixed"
+    meshTypes = MESHTYPES
+    meshDimension = 2
+    expectedOrders = 2
+    meshSizes = [4, 8, 16, 32]
+
+    def getError(self):
+        #Test function
+        phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
+        j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
+        j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
+        phideriv_funX = lambda x: -j_funX(x)
+        phideriv_funY = lambda x: -j_funY(x)
+        q_fun = lambda x: +2*(np.pi**2)*phi_fun(x)
+
+        xc_ana = phi_fun(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        jX_ana = j_funX(self.M.gridFx)
+        jY_ana = j_funY(self.M.gridFy)
+        j_ana = np.r_[jX_ana,jY_ana]
+
+        # Get boundary locations
+        fxm,fxp,fym,fyp = self.M.faceBoundaryInd
+        gBFxm = self.M.gridFx[fxm,:]
+        gBFxp = self.M.gridFx[fxp,:]
+        gBFym = self.M.gridFy[fym,:]
+        gBFyp = self.M.gridFy[fyp,:]
+
+        # Setup Mixed B.C (alpha, beta, gamma)
+        alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        beta_xm, beta_xp =  np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+        beta_ym, beta_yp =  np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+
+        phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
+        phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
+
+        phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
+        phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
+
+        gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
+        gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
+        gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
+        gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
+        gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp)
+
+        alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
+        beta = [beta_xm, beta_xp, beta_ym, beta_yp]
+        gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
+
+        x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
+
+
+        sigma = np.ones(self.M.nC)
+        Mfrho = self.M.getFaceInnerProduct(1./sigma)
+        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
+        V = Utils.sdiag(self.M.vol)
+        Div = V*self.M.faceDiv
+        P_BC, B = self.M.getBCProjWF_simple()
+        q = q_fun(self.M.gridCC)
+        M = B*self.M.aveCC2F
+        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
+        rhs = V*q + Div*MfrhoI*P_BC*x_BC
+        A = Div*MfrhoI*G
+
+        if self.myTest == 'xc':
+            Ainv = Solver(A)
+            xc = Ainv*rhs
+            err = np.linalg.norm((xc-xc_ana), np.inf)
+        else:
+            NotImplementedError
+        return err
+
+
+    def test_order(self):
+        print "==== Testing Mixed boudary conduction for CC-problem ===="
+        self.name = "2D"
+        self.myTest = 'xc'
+        self.orderTest()
+
+class Test3D_InhomogeneousMixed(Tests.OrderTest):
+    name = "3D - Mixed"
+    meshTypes = MESHTYPES
+    meshDimension = 3
+    expectedOrders = 2
+    meshSizes = [4, 8, 16]
+
+    def getError(self):
+        #Test function
+        phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
+        j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
+        j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
+        j_funZ = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.sin(np.pi*x[:,2])
+
+        phideriv_funX = lambda x: -j_funX(x)
+        phideriv_funY = lambda x: -j_funY(x)
+        phideriv_funZ = lambda x: -j_funZ(x)
+
+        q_fun = lambda x: 3*(np.pi**2)*phi_fun(x)
+
+        xc_ana = phi_fun(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        jX_ana = j_funX(self.M.gridFx)
+        jY_ana = j_funY(self.M.gridFy)
+        j_ana = np.r_[jX_ana,jY_ana,jY_ana]
+
+        # Get boundary locations
+        fxm,fxp,fym,fyp,fzm,fzp = self.M.faceBoundaryInd
+        gBFxm = self.M.gridFx[fxm,:]
+        gBFxp = self.M.gridFx[fxp,:]
+        gBFym = self.M.gridFy[fym,:]
+        gBFyp = self.M.gridFy[fyp,:]
+        gBFzm = self.M.gridFz[fzm,:]
+        gBFzp = self.M.gridFz[fzp,:]
+
+        # Setup Mixed B.C (alpha, beta, gamma)
+        alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        beta_xm, beta_xp =  np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+        beta_ym, beta_yp =  np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+        alpha_zm, alpha_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
+        beta_zm, beta_zp =  np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
+
+
+        phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
+        phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
+        phi_bc_zm, phi_bc_zp = phi_fun(gBFzm), phi_fun(gBFzp)
+
+        phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
+        phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
+        phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funZ(gBFzm), phideriv_funZ(gBFzp)
+
+        gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
+        gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
+        gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
+        gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
+        gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp)
+        gamma_zm = gamma_fun(alpha_zm, beta_zm, phi_bc_zm, phiderivY_bc_zm)
+        gamma_zp = gamma_fun(alpha_zp, beta_zp, phi_bc_zp, phiderivY_bc_zp)
+
+        alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
+        beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
+        gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
+
+        x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
+
+
+        sigma = np.ones(self.M.nC)
+        Mfrho = self.M.getFaceInnerProduct(1./sigma)
+        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
+        V = Utils.sdiag(self.M.vol)
+        Div = V*self.M.faceDiv
+        P_BC, B = self.M.getBCProjWF_simple()
+        q = q_fun(self.M.gridCC)
+        M = B*self.M.aveCC2F
+        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
+        rhs = V*q + Div*MfrhoI*P_BC*x_BC
+        A = Div*MfrhoI*G
+
+        if self.myTest == 'xc':
+            #TODO: fix the null space
+            Ainv = Solver(A)
+            xc = Ainv*rhs
+            err = np.linalg.norm((xc-xc_ana), np.inf)
+        else:
+            NotImplementedError
+        return err
+
+
+    def test_order(self):
+        print "==== Testing Mixed boudary conduction for CC-problem ===="
+        self.name = "3D"
+        self.myTest = 'xc'
+        self.orderTest()
+
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/mesh/test_cylMesh.py

@@ -146,6 +146,20 @@ class TestCyl2DMesh(unittest.TestCase):
 
         assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
 
+    def test_getInterpMatCartMesh_Cells2Nodes(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        mc = np.arange(Mc.nC)
+        xr = np.linspace(0,0.4,50)
+        xc = np.linspace(0,0.4,50) + 0.2
+        Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'N')
+        Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC')
+        Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC', locTypeTo='N')
+
+        assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
+
     def test_getInterpMatCartMesh_Faces(self):
 
         Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
@@ -177,6 +191,37 @@ class TestCyl2DMesh(unittest.TestCase):
         assert np.abs(mag[dist > 0.1].min() - 1) < TOL
 
 
+    def test_getInterpMatCartMesh_Faces2Edges(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pf2e = Mc.getInterpolationMatCartMesh(Mr, 'F', locTypeTo='E')
+        mf = np.ones(Mc.nF)
+
+        ecart = Pf2e * mf
+
+        excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey')
+        ezcc = Mr.r(ecart, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 1) < TOL
+        assert np.abs(float(excc[indY]) - 0) < TOL
+        assert np.abs(float(eycc[indX]) - 0) < TOL
+        assert np.abs(float(eycc[indY]) - 1) < TOL
+        assert np.abs((ezcc - 1).sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
     def test_getInterpMatCartMesh_Edges(self):
 
         Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
@@ -185,11 +230,42 @@ class TestCyl2DMesh(unittest.TestCase):
         Pe = Mc.getInterpolationMatCartMesh(Mr, 'E')
         me = np.ones(Mc.nE)
 
-        erect = Pe * me
+        ecart = Pe * me
 
-        excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex')
-        eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey')
-        ezcc = Mr.r(erect, 'E', 'Ez')
+        excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey')
+        ezcc = Mr.aveEz2CC*Mr.r(ecart, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 0) < TOL
+        assert np.abs(float(excc[indY]) + 1) < TOL
+        assert np.abs(float(eycc[indX]) - 1) < TOL
+        assert np.abs(float(eycc[indY]) - 0) < TOL
+        assert np.abs(ezcc.sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
+    def test_getInterpMatCartMesh_Edges2Faces(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pe2f = Mc.getInterpolationMatCartMesh(Mr, 'E', locTypeTo='F')
+        me = np.ones(Mc.nE)
+
+        frect = Pe2f * me
+
+        excc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx')
+        eycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy')
+        ezcc = Mr.r(frect, 'F', 'Fz')
 
         indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
         indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])