Initial merge and minor refactor of simpegPF.

SimPEG/Examples/Mesh_Basic_Types.py

@@ -1,4 +1,5 @@
-from SimPEG import *
+from SimPEG import Mesh, Utils, np
+
 
 def run(plotIt=True):
     """
@@ -8,15 +9,15 @@ def run(plotIt=True):
         Here we show SimPEG used to create three different types of meshes.
 
     """
-    sz = [16,16]
+    sz = [16, 16]
     tM = Mesh.TensorMesh(sz)
     qM = Mesh.TreeMesh(sz)
     qM.refine(lambda cell: 4 if np.sqrt(((np.r_[cell.center]-0.5)**2).sum()) < 0.4 else 3)
-    rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
+    rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz, 'rotate'))
 
     if plotIt:
         import matplotlib.pyplot as plt
-        fig, axes = plt.subplots(1,3,figsize=(14,5))
+        fig, axes = plt.subplots(1, 3, figsize=(14, 5))
         opts = {}
         tM.plotGrid(ax=axes[0], **opts)
         axes[0].set_title('TensorMesh')

SimPEG/Examples/PF_Magnetics_Analytics.py

@@ -0,0 +1,65 @@
+from SimPEG import Mesh, np, PF
+
+
+def run(plotIt=True):
+    """
+        PF: Magnetics: Analytics
+        ========================
+
+        Comparing the magnetics field in Vancouver to Seoul
+
+    """
+
+    xr = np.linspace(-300, 300, 41)
+    yr = np.linspace(-300, 300, 41)
+    X, Y = np.meshgrid(xr, yr)
+    Z = np.ones((np.size(xr), np.size(yr)))*150
+
+    # Bz component in Korea
+    inckr = -8. + 3./60
+    deckr = 54. + 9./60
+    btotkr = 50898.6
+    Bokr = PF.MagAnalytics.IDTtoxyz(inckr, deckr, btotkr)
+
+    bx, by, bz = PF.MagAnalytics.MagSphereAnaFunA(
+        X, Y, Z, 100., 0., 0., 0., 0.01, Bokr, 'secondary'
+    )
+    Bzkr = np.reshape(bz, (np.size(xr), np.size(yr)), order='F')
+
+    # Bz component in Canada
+    incca = 16. + 49./60
+    decca = 70. + 19./60
+    btotca = 54692.1
+    Boca = PF.MagAnalytics.IDTtoxyz(incca, decca, btotca)
+
+    bx, by, bz = PF.MagAnalytics.MagSphereAnaFunA(
+        X, Y, Z, 100., 0., 0., 0., 0.01, Boca, 'secondary'
+    )
+    Bzca = np.reshape(bz, (np.size(xr), np.size(yr)), order='F')
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        from mpl_toolkits.axes_grid1 import make_axes_locatable
+        fig = plt.figure(figsize=(14, 5))
+
+        ax1 = plt.subplot(121)
+        dat1 = plt.imshow(Bzkr, extent=[min(xr), max(xr), min(yr), max(yr)])
+        divider = make_axes_locatable(ax1)
+        cax1 = divider.append_axes("right", size="5%", pad=0.05)
+        ax1.set_xlabel('East-West (m)')
+        ax1.set_ylabel('South-North (m)')
+        plt.colorbar(dat1, cax=cax1)
+        ax1.set_title('$B_z$ field at Seoul, South Korea')
+
+        ax2 = plt.subplot(122)
+        dat2 = plt.imshow(Bzca, extent=[min(xr), max(xr), min(yr), max(yr)])
+        divider = make_axes_locatable(ax2)
+        cax2 = divider.append_axes("right", size="5%", pad=0.05)
+        ax2.set_xlabel('East-West (m)')
+        ax2.set_ylabel('South-North (m)')
+        plt.colorbar(dat2, cax=cax2)
+        ax2.set_title('$B_z$ field at Vancouver, Canada')
+        plt.show()
+
+if __name__ == '__main__':
+    run()

SimPEG/Examples/__init__.py

@@ -20,9 +20,10 @@ import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
 import MT_1D_ForwardAndInversion
 import MT_3D_Foward
+import PF_Magnetics_Analytics
 import Utils_surface2ind_topo
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "PF_Magnetics_Analytics", "Utils_surface2ind_topo"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/PF/BaseGrav.py

@@ -0,0 +1,53 @@
+from SimPEG import Maps, Survey, Utils, np, sp
+from scipy.constants import mu_0
+import re
+
+
+class LinearSurvey(Survey.BaseSurvey):
+    """Base Magnetics Survey"""
+
+    rxLoc = None #: receiver locations
+    rxType = None #: receiver type
+
+    def __init__(self, srcField, **kwargs):
+        self.srcField = srcField
+        Survey.BaseSurvey.__init__(self, **kwargs)
+
+    def eval(self, u):
+        return u
+
+    @property
+    def nD(self):
+        return self.prob.G.shape[0]
+
+    @property
+    def nRx(self):
+        return self.srcField.rxList[0].locs.shape[0]
+    # def setBackgroundField(self, SrcField):
+
+    #     if getattr(self, 'B0', None) is None:
+    #         self._B0 = SrcField.param[0] * dipazm_2_xyz(SrcField.param[1],SrcField.param[2])
+
+    #     return self._B0
+
+
+class SrcField(Survey.BaseSrc):
+    """ Define the inducing field """
+
+    param = None #: Inducing field param (Amp, Incl, Decl)
+
+    def __init__(self, rxList, **kwargs):
+        super(SrcField, self).__init__(rxList, **kwargs)
+
+
+class RxObs(Survey.BaseRx):
+    """A station location must have be located in 3-D"""
+    def __init__(self, locsXYZ, **kwargs):
+        locs = locsXYZ
+        assert locsXYZ.shape[1] == 3, 'locs must in 3-D (x,y,z).'
+        super(RxObs, self).__init__(locs, 'tmi', storeProjections=False, **kwargs)
+
+    @property
+    def nD(self):
+        """Number of data in the receiver."""
+        return self.locs[0].shape[0]

SimPEG/PF/BaseMag.py

@@ -0,0 +1,194 @@
+from SimPEG import Maps, Survey, Utils, np, sp
+from scipy.constants import mu_0
+import re
+
+
+class BaseMagSurvey(Survey.BaseSurvey):
+    """Base Magnetics Survey"""
+
+    rxLoc = None   #: receiver locations
+    rxType = None   #: receiver type
+
+    def __init__(self, **kwargs):
+        Survey.BaseSurvey.__init__(self, **kwargs)
+
+    def setBackgroundField(self, Inc, Dec, Btot):
+
+        Bx = Btot*np.cos(Inc/180.*np.pi)*np.sin(Dec/180.*np.pi)
+        By = Btot*np.cos(Inc/180.*np.pi)*np.cos(Dec/180.*np.pi)
+        Bz = -Btot*np.sin(Inc/180.*np.pi)
+
+        self.B0 = np.r_[Bx, By, Bz]
+
+    @property
+    def Qfx(self):
+        if getattr(self, '_Qfx', None) is None:
+            self._Qfx = self.prob.mesh.getInterpolationMat(self.rxLoc, 'Fx')
+        return self._Qfx
+
+    @property
+    def Qfy(self):
+        if getattr(self, '_Qfy', None) is None:
+            self._Qfy = self.prob.mesh.getInterpolationMat(self.rxLoc, 'Fy')
+        return self._Qfy
+
+    @property
+    def Qfz(self):
+        if getattr(self, '_Qfz', None) is None:
+            self._Qfz = self.prob.mesh.getInterpolationMat(self.rxLoc, 'Fz')
+        return self._Qfz
+
+    def projectFields(self, u):
+        """
+            This function projects the fields onto the data space.
+
+            Especially, here for we use total magnetic intensity (TMI) data,
+            which is common in practice.
+
+            First we project our B on to data location
+
+            .. math::
+
+                \mathbf{B}_{rec} = \mathbf{P} \mathbf{B}
+
+            then we take the dot product between B and b_0
+
+            .. math ::
+
+                \\text{TMI} = \\vec{B}_s \cdot \hat{B}_0
+
+        """
+        # TODO: There can be some different tyes of data like |B| or B
+
+        bfx = self.Qfx*u['B']
+        bfy = self.Qfy*u['B']
+        bfz = self.Qfz*u['B']
+
+        # Generate unit vector
+        B0 = self.prob.survey.B0
+        Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2)
+        box = B0[0]/Bot
+        boy = B0[1]/Bot
+        boz = B0[2]/Bot
+
+        # return bfx*box + bfx*boy + bfx*boz
+        return bfx*box + bfy*boy + bfz*boz
+
+    @Utils.count
+    def projectFieldsDeriv(self, B):
+        """
+            This function projects the fields onto the data space.
+
+            .. math::
+
+                \\frac{\partial d_\\text{pred}}{\partial \mathbf{B}} = \mathbf{P}
+
+            Especially, this function is for TMI data type
+
+        """
+        # Generate unit vector
+        B0 = self.prob.survey.B0
+        Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2)
+        box = B0[0]/Bot
+        boy = B0[1]/Bot
+        boz = B0[2]/Bot
+
+        return self.Qfx*box+self.Qfy*boy+self.Qfz*boz
+
+    def projectFieldsAsVector(self, B):
+
+        bfx = self.Qfx*B
+        bfy = self.Qfy*B
+        bfz = self.Qfz*B
+
+        return np.r_[bfx, bfy, bfz]
+
+
+class LinearSurvey(Survey.BaseSurvey):
+    """Base Magnetics Survey"""
+
+    rxLoc = None #: receiver locations
+    rxType = None #: receiver type
+
+    def __init__(self, srcField, **kwargs):
+        self.srcField = srcField
+        Survey.BaseSurvey.__init__(self, **kwargs)
+
+    def eval(self, u):
+        return u
+
+    @property
+    def nD(self):
+        return self.prob.G.shape[0]
+
+    @property
+    def nRx(self):
+        return self.srcField.rxList[0].locs.shape[0]
+    # def setBackgroundField(self, SrcField):
+
+    #     if getattr(self, 'B0', None) is None:
+    #         self._B0 = SrcField.param[0] * dipazm_2_xyz(SrcField.param[1],SrcField.param[2])
+
+    #     return self._B0
+
+
+class SrcField(Survey.BaseSrc):
+    """ Define the inducing field """
+
+    param = None #: Inducing field param (Amp, Incl, Decl)
+
+    def __init__(self, rxList, **kwargs):
+        super(SrcField, self).__init__(rxList, **kwargs)
+
+
+class RxObs(Survey.BaseRx):
+    """A station location must have be located in 3-D"""
+    def __init__(self, locsXYZ, **kwargs):
+        locs = locsXYZ
+        assert locsXYZ.shape[1] == 3, 'locs must in 3-D (x,y,z).'
+        super(RxObs, self).__init__(locs, 'tmi', storeProjections=False, **kwargs)
+
+    @property
+    def nD(self):
+        """Number of data in the receiver."""
+        return self.locs[0].shape[0]
+
+
+class MagSurveyBx(object):
+    """docstring for MagSurveyBx"""
+    def __init__(self, **kwargs):
+        Survey.BaseData.__init__(self, **kwargs)
+
+    def projectFields(self, B):
+        bfx = self.Qfx*B
+        return bfx
+
+
+class BaseMagMap(Maps.IdentityMap):
+    """BaseMagMap"""
+
+    def __init__(self, mesh, **kwargs):
+        Maps.IdentityMap.__init__(self, mesh)
+
+    def _transform(self, m):
+
+        return mu_0*(1 + m)
+
+    def deriv(self, m):
+
+        return mu_0*sp.identity(self.nP)
+
+
+class WeightMap(Maps.IdentityMap):
+    """Weighted Map for distributed parameters"""
+
+    def __init__(self, nP, weight, **kwargs):
+        Maps.IdentityMap.__init__(self, nP)
+        self.mesh = None
+        self.weight = weight
+
+    def _transform(self, m):
+        return m*self.weight
+
+    def deriv(self, m):
+        return Utils.sdiag(self.weight)

SimPEG/PF/Gravity.py

@@ -0,0 +1,491 @@
+from SimPEG import *
+import BaseGrav as GRAV
+import re
+
+
+class GravityIntegral(Problem.BaseProblem):
+
+    # surveyPair = Survey.LinearSurvey
+
+    storeG = True #: Store the forward matrix by default, otherwise just compute d
+    actInd = None #: Active cell indices provided
+
+    def __init__(self, mesh, mapping=None, **kwargs):
+        Problem.BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs)
+
+    def fwr_op(self):
+        # Add forward function
+        # kappa = self.curModel.kappa TODO
+        sus = self.mapping*self.curModel
+        return self.G.dot(sus)
+
+    def fields(self, m):
+        self.curModel = m
+        total = np.zeros(self.survey.nRx)
+        induced = self.fwr_op()
+        # rem = self.rem
+
+        if induced is not None:
+            total += induced
+
+        return total
+
+        # return self.G.dot(self.mapping*(m))
+
+    def Jvec(self, m, v, f=None):
+        dmudm = self.mapping.deriv(m)
+        return self.G.dot(dmudm*v)
+
+    def Jtvec(self, m, v, f=None):
+        dmudm = self.mapping.deriv(m)
+        return dmudm.T * (self.G.T.dot(v))
+
+    @property
+    def G(self):
+        if not self.ispaired:
+            raise Exception('Need to pair!')
+
+        if getattr(self, '_G', None) is None:
+            self._G = self.Intrgl_Fwr_Op( 'z' )
+
+        return self._G
+
+    def Intrgl_Fwr_Op(self, flag):
+
+        """
+
+        Gravity forward operator in integral form
+
+        flag        = 'z' | 'xyz'
+
+        Return
+        _G        = Linear forward modeling operation
+
+        Created on March, 15th 2016
+
+        @author: dominiquef
+
+         """
+        # Find non-zero cells
+        # inds = np.nonzero(actv)[0]
+        if getattr(self, 'actInd', None) is not None:
+
+            if self.actInd.dtype=='bool':
+                inds = np.asarray([inds for inds, elem in enumerate(self.actInd, 1) if elem], dtype = int) - 1
+            else:
+                inds = self.actInd
+
+        else:
+
+            inds = np.asarray(range(self.mesh.nC))
+
+        nC = len(inds)
+
+        # Create active cell projector
+        P = sp.csr_matrix(
+            (np.ones(nC), (inds, range(nC))),
+            shape=(self.mesh.nC, nC)
+        )
+
+        # Create vectors of nodal location (lower and upper corners for each cell)
+        xn = self.mesh.vectorNx
+        yn = self.mesh.vectorNy
+        zn = self.mesh.vectorNz
+
+        yn2, xn2, zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:])
+        yn1, xn1, zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1])
+
+        Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)]
+        Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)]
+        Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)]
+
+        rxLoc = self.survey.srcField.rxList[0].locs
+        ndata = rxLoc.shape[0]
+
+        # Pre-allocate space and create magnetization matrix if required
+        # Pre-allocate space
+        if flag == 'z':
+
+            G = np.zeros((ndata, nC))
+
+        elif flag == 'xyz':
+
+            G = np.zeros((int(3*ndata), nC))
+
+        else:
+
+            print """Flag must be either 'z' | 'xyz', please revised"""
+            return
+
+
+        # Loop through all observations and create forward operator (ndata-by-nC)
+        print "Begin calculation of forward operator: " + flag
+
+        # Add counter to dsiplay progress. Good for large problems
+        count = -1;
+        for ii in range(ndata):
+
+            if flag=='z':
+                tt = get_T_mat(Xn, Yn, Zn, rxLoc[ii, :])
+                G[ii, :] = tt
+
+            elif flag == 'xyz':
+                print "Sorry 3-component not implemented yet"
+
+            # Display progress
+            count = progress(ii, count, ndata)
+
+        print "Done 100% ...forward operator completed!!\n"
+
+        return G
+
+
+def get_T_mat(Xn, Yn, Zn, rxLoc):
+    """
+    Load in the active nodes of a tensor mesh and computes the gravity tensor
+    for a given observation location rxLoc[obsx, obsy, obsz]
+
+    INPUT:
+    Xn, Yn, Zn: Node location matrix for the lower and upper most corners of
+                all cells in the mesh shape[nC,2]
+    M
+    OUTPUT:
+    Tx = [Txx Txy Txz]
+    Ty = [Tyx Tyy Tyz]
+    Tz = [Tzx Tzy Tzz]
+
+    where each elements have dimension 1-by-nC.
+    Only the upper half 5 elements have to be computed since symetric.
+    Currently done as for-loops but will eventually be changed to vector
+    indexing, once the topography has been figured out.
+
+    """
+    NewtG=6.6738e-3
+    eps = 1e-10 # add a small value to the locations to avoid /0
+
+    nC = Xn.shape[0]
+
+    # Pre-allocate space for 1D array
+    T = np.zeros((1,nC))
+
+    dz = rxLoc[2] - Zn + eps
+
+    dy = Yn - rxLoc[1] + eps
+
+    dx = Xn - rxLoc[0] + eps
+
+    # Compute contribution from each corners
+    for aa in range(2):
+        for bb in range(2):
+            for cc in range(2):
+
+                r = (
+                        dx[:, aa] ** 2 +
+                        dy[:, bb] ** 2 +
+                        dz[:, cc] ** 2
+                    ) ** (0.50)
+
+                T = T - NewtG * (-1) ** aa * (-1) ** bb * (-1) ** cc * (
+                    dx[:, aa] * np.log(dy[:, bb] + r) +
+                    dy[:, bb] * np.log(dx[:, aa] + r) -
+                    dz[:, cc] * np.arctan(
+                        dx[:, aa] * dy[:, bb] / (dz[:, cc] * r)
+                    )
+                )
+
+    return T
+
+
+def progress(iter, prog, final):
+    """
+    progress(iter,prog,final)
+
+    Function measuring the progress of a process and print to screen the %.
+    Useful to estimate the remaining runtime of a large problem.
+
+    Created on Dec, 20th 2015
+
+    @author: dominiquef
+    """
+    arg = np.floor(float(iter)/float(final)*10.)
+
+    if  arg > prog:
+
+        strg = "Done " + str(arg*10) + " %"
+        print strg
+        prog = arg
+
+    return prog
+
+
+def writeUBCobs(filename, survey, d):
+    """
+    writeUBCobs(filename,survey,d)
+
+    Function writing an observation file in UBC-GRAV3D format.
+
+    INPUT
+    filename    : Name of out file including directory
+    survey
+    flag          : dobs | dpred
+
+    OUTPUT
+    Obsfile
+
+    """
+
+    rxLoc = survey.srcField.rxList[0].locs
+
+    wd = survey.std
+
+    data = np.c_[rxLoc , d , wd]
+
+    with file(filename,'w') as fid:
+        fid.write('%i\n' %len(d) )
+        np.savetxt(fid, data, fmt='%e', delimiter=' ', newline='\n')
+
+
+    print "Observation file saved to: " + filename
+
+
+def getActiveTopo(mesh, topo, flag):
+    """
+    getActiveTopo(mesh,topo)
+
+    Function creates an active cell model from topography
+
+    INPUT
+    mesh        : Mesh in SimPEG format
+    topo        : Scatter points defining topography [x,y,z]
+
+    OUTPUT
+    actv        : Active cell model
+
+    """
+    import scipy.interpolate as interpolation
+
+    if flag == 'N':
+        Zn = np.zeros((mesh.nNx, mesh.nNy))
+        # wght    = np.zeros((mesh.nNx,mesh.nNy))
+        cx = mesh.vectorNx
+        cy = mesh.vectorNy
+
+    F = interpolation.NearestNDInterpolator(topo[:, 0:2], topo[:, 2])
+    [Y, X] = np.meshgrid(cy, cx)
+
+    Zn = F(X, Y)
+
+    actv = np.zeros((mesh.nCx, mesh.nCy, mesh.nCz))
+
+    if flag == 'N':
+        Nz = mesh.vectorNz[1:]
+
+    for jj in range(mesh.nCy):
+
+        for ii in range(mesh.nCx):
+
+            temp = [kk for kk in range(len(Nz)) if np.all(Zn[ii:(ii+2), jj:(jj+2)] > Nz[kk]) ]
+            actv[ii, jj, temp] = 1
+
+    actv = mkvc(actv == 1)
+
+    inds = np.asarray([inds for inds, elem in enumerate(actv, 1) if elem], dtype = int) - 1
+
+    return inds
+
+def plot_obs_2D(survey,varstr):
+    """ Function plot_obs(rxLoc,d,wd)
+    Generate a 2d interpolated plot from scatter points of data
+
+    INPUT
+    rxLoc       : Observation locations [x,y,z]
+    d           : Data vector
+    wd          : Uncertainty vector
+
+    OUTPUT
+    figure()
+
+    Created on Dec, 27th 2015
+
+    @author: dominiquef
+
+    """
+
+    from scipy.interpolate import griddata
+    import pylab as plt
+
+    rxLoc   = survey.srcField.rxList[0].locs
+    d       = survey.dobs
+    wd      = survey.std
+
+    # Create grid of points
+    x = np.linspace(rxLoc[:,0].min(), rxLoc[:,0].max(), 100)
+    y = np.linspace(rxLoc[:,1].min(), rxLoc[:,1].max(), 100)
+
+    X, Y = np.meshgrid(x,y)
+
+    # Interpolate
+    d_grid = griddata(rxLoc[:,0:2],d,(X,Y), method ='linear')
+
+    # Plot result
+    plt.figure()
+    plt.subplot()
+    plt.imshow(d_grid, extent=[x.min(), x.max(), y.min(), y.max()],origin = 'lower')
+    plt.colorbar(fraction=0.02)
+    plt.contour(X,Y, d_grid,10)
+    plt.scatter(rxLoc[:,0],rxLoc[:,1], c=d, s=20)
+    plt.title(varstr)
+    plt.gca().set_aspect('equal', adjustable='box')
+
+def readUBCgravObs(obs_file):
+
+    """
+    Read UBC grav file format
+
+    INPUT:
+    :param fileName, path to the UBC obs grav file
+
+    OUTPUT:
+    :param survey
+
+    """
+
+    fid = open(obs_file,'r')
+
+    # First line has the number of rows
+    line = fid.readline()
+    ndat = np.array(line.split(),dtype=int)
+
+    # Pre-allocate space for obsx, obsy, obsz, data, uncert
+    line = fid.readline()
+    temp = np.array(line.split(),dtype=float)
+
+    d  = np.zeros(ndat, dtype=float)
+    wd = np.zeros(ndat, dtype=float)
+    locXYZ = np.zeros( (ndat,3), dtype=float)
+
+    for ii in range(ndat):
+
+        temp = np.array(line.split(),dtype=float)
+        locXYZ[ii,:] = temp[:3]
+        d[ii] = temp[3]
+        wd[ii] = temp[4]
+        line = fid.readline()
+
+    rxLoc = GRAV.RxObs(locXYZ)
+    srcField = GRAV.SrcField([rxLoc])
+    survey = GRAV.LinearSurvey(srcField)
+    survey.dobs =  d
+    survey.std =  wd
+    return survey
+
+
+def read_GRAVinv_inp(input_file):
+    """Read input files for forward modeling MAG data with integral form
+    INPUT:
+    input_file: File name containing the forward parameter
+
+    OUTPUT:
+    mshfile
+    obsfile
+    topofile
+    start model
+    ref model
+    weightfile
+    chi_target
+    as, ax ,ay, az
+    upper, lower bounds
+    lp, lqx, lqy, lqz
+
+    # All files should be in the working directory, otherwise the path must
+    # be specified.
+
+    Created on Dec 21th, 2015
+
+    @author: dominiquef
+    """
+
+
+    fid = open(input_file,'r')
+
+    # Line 1
+    line = fid.readline()
+    l_input  = line.split('!')
+    mshfile = l_input[0].rstrip()
+
+    # Line 2
+    line = fid.readline()
+    l_input  = line.split('!')
+    obsfile = l_input[0].rstrip()
+
+    # Line 3
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input=='null':
+        topofile = []
+
+    else:
+        topofile = l_input[0].rstrip()
+
+
+    # Line 4
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input[0]=='VALUE':
+        mstart = float(l_input[1])
+
+    else:
+        mstart = l_input[0].rstrip()
+
+    # Line 5
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input[0]=='VALUE':
+        mref = float(l_input[1])
+
+    else:
+        mref = l_input[0].rstrip()
+
+
+    # Line 7
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input[0]=='DEFAULT':
+        wgtfile = None
+
+    else:
+        wgtfile = l_input[0].rstrip()
+
+    # Line 8
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    chi = float(l_input[0])
+
+    # Line 9
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    val = np.array(l_input[0:4])
+    alphas = val.astype(np.float)
+
+    # Line 10
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input[0]=='VALUE':
+        val   = np.array(l_input[1:3])
+        bounds = val.astype(np.float)
+
+    else:
+        bounds = l_input[0].rstrip()
+
+    # Line 11
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input[0]=='VALUE':
+        val   = np.array(l_input[1:6])
+        lpnorms = val.astype(np.float)
+
+    else:
+        lpnorms = l_input[0].rstrip()
+
+    return mshfile, obsfile, topofile, mstart, mref, wgtfile, chi, alphas, bounds, lpnorms
+

SimPEG/PF/GravityDriver.py

@@ -0,0 +1,295 @@
+import re, os
+from SimPEG import Mesh, np, Utils
+import BaseGrav, Gravity
+
+
+class GravityDriver_Inv(object):
+    """docstring for GravityDriver_Inv"""
+
+    def __init__(self, input_file=None):
+        if input_file is not None:
+            self.basePath = os.path.sep.join(input_file.split(os.path.sep)[:-1])
+            if len(self.basePath) > 0:
+                self.basePath += os.path.sep
+            self.readDriverFile(input_file.split(os.path.sep)[-1])
+
+    def readDriverFile(self, input_file):
+        """
+            Read input files for forward modeling GRAV data with integral form
+            INPUT:
+            input_file: File name containing the forward parameter
+
+            OUTPUT:
+            mshfile
+            obsfile
+            topofile
+            start model
+            ref model
+            active cells model
+            weightfile
+            chi_target
+            as, ax ,ay, az
+            upper, lower bounds
+            lp, lqx, lqy, lqz
+            eps_p, eps_q
+            # All files should be in the working directory, otherwise the path must
+            # be specified.
+
+        """
+
+        fid = open(self.basePath + input_file, 'r')
+
+        # Line 1
+        line = fid.readline()
+        l_input  = line.split('!')
+        mshfile = l_input[0].rstrip()
+
+        # Line 2
+        line = fid.readline()
+        l_input  = line.split('!')
+        obsfile = l_input[0].rstrip()
+
+        # Line 3
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input=='null':
+            topofile = []
+
+        else:
+            topofile = l_input[0].rstrip()
+
+
+        # Line 4
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input[0]=='VALUE':
+            mstart = float(l_input[1])
+
+        else:
+            mstart = l_input[0].rstrip()
+
+        # Line 5
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            mref = float(l_input[1])
+
+        else:
+            mref = l_input[0].rstrip()
+
+        # Line 6
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input[0]=='VALUE':
+            staticInput = float(l_input[1])
+
+        elif l_input[0]=='DEFAULT':
+            staticInput = None
+
+        else:
+            staticInput = l_input[0].rstrip()
+
+        # Line 7
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input=='DEFAULT':
+            wgtfile = []
+
+        else:
+            wgtfile = l_input[0].rstrip()
+
+        # Line 8
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        chi = float(l_input[0])
+
+        # Line 9
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        val = np.array(l_input[0:4])
+        alphas = val.astype(np.float)
+
+        # Line 10
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:3])
+            bounds = val.astype(np.float)
+
+        else:
+            bounds = l_input[0].rstrip()
+
+        # Line 11
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:6])
+            lpnorms = val.astype(np.float)
+
+        else:
+            lpnorms = l_input[0].rstrip()
+
+        # Line 12
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:3])
+            eps = val.astype(np.float)
+
+        else:
+            eps = [None, None]
+
+        self.mshfile  = mshfile
+        self.obsfile  = obsfile
+        self.topofile = topofile
+        self.mstart   = mstart
+        self._mrefInput = mref
+        self._staticInput = staticInput
+        self.wgtfile  = wgtfile
+        self.chi      = chi
+        self.alphas   = alphas
+        self.bounds   = bounds
+        self.lpnorms  = lpnorms
+        self.eps      = eps
+
+    @property
+    def mesh(self):
+        if getattr(self, '_mesh', None) is None:
+            self._mesh = Mesh.TensorMesh.readUBC(self.basePath + self.mshfile)
+        return self._mesh
+
+    @property
+    def survey(self):
+        if getattr(self, '_survey', None) is None:
+            self._survey = self.readGravityObservations(self.basePath + self.obsfile)
+        return self._survey
+
+    @property
+    def activeCells(self):
+        if getattr(self, '_activeCells', None) is None:
+            if self.topofile == 'null':
+                self._activeCells = np.arange(mesh.nC)
+            else:
+                topo = np.genfromtxt(self.basePath + self.topofile, skip_header=1)
+                # Find the active cells
+                active = Utils.surface2ind_topo(self.mesh,topo,'N')
+                inds = np.asarray([inds for inds, elem in enumerate(active, 1) if elem], dtype = int) - 1
+                self._activeCells = inds
+
+        return self._activeCells
+
+    @property
+    def staticCells(self):
+        if getattr(self, '_staticCells', None) is None:
+
+            if getattr(self, '_staticInput', None) is None:
+                # All cells are dynamic: 1's
+                self._dynamicCells = np.arange(len(self.m0))
+                self._staticCells = []
+
+            # Cells with specific value are static: 0's
+            else:
+                if isinstance(self._staticInput, float):
+                    staticCells = self.m0 == self._staticInput
+
+                else:
+                    # Read from file active cells with 0:air, 1:dynamic, -1 static
+                    staticCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
+                    staticCells = staticCells[self.activeCells] == -1
+
+                inds = np.asarray([inds for inds, elem in enumerate(staticCells, 1) if elem], dtype = int) - 1
+                self._staticCells = inds
+
+        return self._staticCells
+
+    @property
+    def dynamicCells(self):
+        if getattr(self, '_dynamicCells', None) is None:
+
+            if getattr(self, '_staticInput', None) is None:
+                # All cells are dynamic: 1's
+                self._dynamicCells = np.arange(len(self.m0))
+
+            # Cells with specific value are static: 0's
+            else:
+                if isinstance(self._staticInput, float):
+                    dynamicCells = self.m0 != self._staticInput
+
+                else:
+                    # Read from file active cells with 0:air, 1:dynamic, -1 static
+                    dynamicCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
+                    dynamicCells = dynamicCells[self.activeCells] == 1
+
+                inds = np.asarray([inds for inds, elem in enumerate(dynamicCells, 1) if elem], dtype = int) - 1
+                self._dynamicCells = inds
+
+        return self._dynamicCells
+
+    @property
+    def nC(self):
+        if getattr(self, '_nC', None) is None:
+            self._nC = len(self.activeCells)
+        return self._nC
+
+    @property
+    def m0(self):
+        if getattr(self, '_m0', None) is None:
+            if isinstance(self.mstart, float):
+                self._m0 = np.ones(self.nC) * self.mstart
+            else:
+
+                self._m0 = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self.mstart)
+                self._m0 = self._m0[self.activeCells]
+
+        return self._m0
+
+    @property
+    def mref(self):
+        if getattr(self, '_mref', None) is None:
+            if isinstance(self._mrefInput, float):
+                self._mref = np.ones(self.nC) * self._mrefInput
+            else:
+                self._mref = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._mrefInput)
+                self._mref = self._mref[self.activeCells]
+        return self._mref
+
+    def readGravityObservations(self, obs_file):
+        """
+        Read UBC grav file format
+
+        INPUT:
+        :param fileName, path to the UBC obs grav file
+
+        OUTPUT:
+        :param survey
+
+        """
+
+        fid = open(obs_file,'r')
+
+        # First line has the number of rows
+        line = fid.readline()
+        ndat = np.array(line.split(),dtype=int)
+
+        # Pre-allocate space for obsx, obsy, obsz, data, uncert
+        line = fid.readline()
+        temp = np.array(line.split(),dtype=float)
+
+        d  = np.zeros(ndat, dtype=float)
+        wd = np.zeros(ndat, dtype=float)
+        locXYZ = np.zeros( (ndat,3), dtype=float)
+
+        for ii in range(ndat):
+
+            temp = np.array(line.split(),dtype=float)
+            locXYZ[ii,:] = temp[:3]
+            d[ii] = temp[3]
+            wd[ii] = temp[4]
+            line = fid.readline()
+
+        rxLoc = BaseGrav.RxObs(locXYZ)
+        srcField = BaseGrav.SrcField([rxLoc])
+        survey = BaseGrav.LinearSurvey(srcField)
+        survey.dobs =  d
+        survey.std =  wd
+        return survey

SimPEG/PF/MagAnalytics.py

@@ -0,0 +1,278 @@
+from scipy.constants import mu_0
+from SimPEG import *
+from SimPEG.Utils import kron3, speye, sdiag
+import matplotlib.pyplot as plt
+
+
+def spheremodel(mesh, x0, y0, z0, r):
+    """
+        Generate model indicies for sphere
+        - (x0, y0, z0 ): is the center location of sphere
+        - r: is the radius of the sphere
+        - it returns logical indicies of cell-center model
+    """
+    ind = np.sqrt( (mesh.gridCC[:,0]-x0)**2+(mesh.gridCC[:,1]-y0)**2+(mesh.gridCC[:,2]-z0)**2 ) < r
+    return ind
+
+
+def MagSphereAnaFun(x, y, z, R, x0, y0, z0, mu1, mu2, H0, flag='total'):
+    """
+        test
+        Analytic function for Magnetics problem. The set up here is
+        magnetic sphere in whole-space assuming that the inducing field is oriented in the x-direction.
+
+        * (x0,y0,z0)
+        * (x0, y0, z0 ): is the center location of sphere
+        * r: is the radius of the sphere
+
+    .. math::
+
+        \mathbf{H}_0 = H_0\hat{x}
+
+
+    """
+
+    if (~np.size(x)==np.size(y)==np.size(z)):
+        print "Specify same size of x, y, z"
+        return
+    dim = x.shape
+    x = Utils.mkvc(x)
+    y = Utils.mkvc(y)
+    z = Utils.mkvc(z)
+
+    ind = np.sqrt((x-x0)**2+(y-y0)**2+(z-z0)**2 ) < R
+    r = Utils.mkvc(np.sqrt((x-x0)**2+(y-y0)**2+(z-z0)**2 ))
+    Bx = np.zeros(x.size)
+    By = np.zeros(x.size)
+    Bz = np.zeros(x.size)
+
+    # Inside of the sphere
+    rf2 = 3*mu1/(mu2+2*mu1)
+    if flag is 'total' and any(ind):
+        Bx[ind] = mu2*H0*(rf2)
+    elif (flag == 'secondary'):
+        Bx[ind] = mu2*H0*(rf2)-mu1*H0
+
+    By[ind] = 0.
+    Bz[ind] = 0.
+    # Outside of the sphere
+    rf1 = (mu2-mu1)/(mu2+2*mu1)
+    if (flag == 'total'):
+        Bx[~ind] = mu1*(H0+H0/r[~ind]**5*(R**3)*rf1*(2*(x[~ind]-x0)**2-(y[~ind]-y0)**2-(z[~ind]-z0)**2))
+    elif (flag == 'secondary'):
+        Bx[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(2*(x[~ind]-x0)**2-(y[~ind]-y0)**2-(z[~ind]-z0)**2))
+
+    By[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*(x[~ind]-x0)*(y[~ind]-y0)))
+    Bz[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*(x[~ind]-x0)*(z[~ind]-z0)))
+    return np.reshape(Bx, x.shape, order='F'), np.reshape(By, x.shape, order='F'), np.reshape(Bz, x.shape, order='F')
+
+
+def CongruousMagBC(mesh, Bo, chi):
+    """
+        Computing boundary condition using Congrous sphere method.
+        This is designed for secondary field formulation.
+
+        >> Input
+
+        * mesh:   Mesh class
+        * Bo:     np.array([Box, Boy, Boz]): Primary magnetic flux
+        * chi:    susceptibility at cell volume
+
+        .. math::
+
+            \\vec{B}(r) = \\frac{\mu_0}{4\pi} \\frac{m}{ \| \\vec{r} - \\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
+
+    """
+
+    ind = chi > 0.
+    V = mesh.vol[ind].sum()
+
+    gamma = 1/V*(chi*mesh.vol).sum() # like a mass!
+
+    Bot = np.sqrt(sum(Bo**2))
+    mx = Bo[0]/Bot
+    my = Bo[1]/Bot
+    mz = Bo[2]/Bot
+
+    mom = 1/mu_0*Bot*gamma*V/(1+gamma/3)
+    xc = sum(chi[ind]*mesh.gridCC[:,0][ind])/sum(chi[ind])
+    yc = sum(chi[ind]*mesh.gridCC[:,1][ind])/sum(chi[ind])
+    zc = sum(chi[ind]*mesh.gridCC[:,2][ind])/sum(chi[ind])
+
+    indxd, indxu, indyd, indyu, indzd, indzu =  mesh.faceBoundaryInd
+
+    const = mu_0/(4*np.pi)*mom
+    rfun = lambda x: np.sqrt((x[:,0]-xc)**2 + (x[:,1]-yc)**2 + (x[:,2]-zc)**2)
+
+    mdotrx = (mx*(mesh.gridFx[(indxd|indxu),0]-xc)/rfun(mesh.gridFx[(indxd|indxu),:]) +
+              my*(mesh.gridFx[(indxd|indxu),1]-yc)/rfun(mesh.gridFx[(indxd|indxu),:]) +
+              mz*(mesh.gridFx[(indxd|indxu),2]-zc)/rfun(mesh.gridFx[(indxd|indxu),:]))
+
+    Bbcx = const/(rfun(mesh.gridFx[(indxd|indxu),:])**3)*(3*mdotrx*(mesh.gridFx[(indxd|indxu),0]-xc)/rfun(mesh.gridFx[(indxd|indxu),:])-mx)
+
+    mdotry = (mx*(mesh.gridFy[(indyd|indyu),0]-xc)/rfun(mesh.gridFy[(indyd|indyu),:]) +
+              my*(mesh.gridFy[(indyd|indyu),1]-yc)/rfun(mesh.gridFy[(indyd|indyu),:]) +
+              mz*(mesh.gridFy[(indyd|indyu),2]-zc)/rfun(mesh.gridFy[(indyd|indyu),:]))
+
+    Bbcy = const/(rfun(mesh.gridFy[(indyd|indyu),:])**3)*(3*mdotry*(mesh.gridFy[(indyd|indyu),1]-yc)/rfun(mesh.gridFy[(indyd|indyu),:])-my)
+
+    mdotrz = (mx*(mesh.gridFz[(indzd|indzu),0]-xc)/rfun(mesh.gridFz[(indzd|indzu),:])  +
+              my*(mesh.gridFz[(indzd|indzu),1]-yc)/rfun(mesh.gridFz[(indzd|indzu),:]) +
+              mz*(mesh.gridFz[(indzd|indzu),2]-zc)/rfun(mesh.gridFz[(indzd|indzu),:]))
+
+    Bbcz = const/(rfun(mesh.gridFz[(indzd|indzu),:])**3)*(3*mdotrz*(mesh.gridFz[(indzd|indzu),2]-zc)/rfun(mesh.gridFz[(indzd|indzu),:])-mz)
+
+    return np.r_[Bbcx, Bbcy, Bbcz], (1/gamma-1/(3+gamma))*1/V
+
+
+def MagSphereAnaFunA(x, y, z, R, xc, yc, zc, chi, Bo, flag):
+    """
+        Computing boundary condition using Congrous sphere method.
+        This is designed for secondary field formulation.
+        >> Input
+        mesh:   Mesh class
+        Bo:     np.array([Box, Boy, Boz]): Primary magnetic flux
+        Chi:    susceptibility at cell volume
+
+        .. math::
+
+            \\vec{B}(r) = \\frac{\mu_0}{4\pi}\\frac{m}{\| \\vec{r}-\\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
+
+    """
+    if (~np.size(x)==np.size(y)==np.size(z)):
+        print "Specify same size of x, y, z"
+        return
+    dim = x.shape
+    x = Utils.mkvc(x)
+    y = Utils.mkvc(y)
+    z = Utils.mkvc(z)
+
+    Bot = np.sqrt(sum(Bo**2))
+    mx = Bo[0]/Bot
+    my = Bo[1]/Bot
+    mz = Bo[2]/Bot
+
+    ind = np.sqrt((x-xc)**2+(y-yc)**2+(z-zc)**2 ) < R
+
+    Bx = np.zeros(x.size)
+    By = np.zeros(x.size)
+    Bz = np.zeros(x.size)
+
+    # Inside of the sphere
+    rf2 = 3/(chi+3)*(1+chi)
+    if (flag == 'total'):
+        Bx[ind] = Bo[0]*(rf2)
+        By[ind] = Bo[1]*(rf2)
+        Bz[ind] = Bo[2]*(rf2)
+    elif (flag == 'secondary'):
+        Bx[ind] = Bo[0]*(rf2)-Bo[0]
+        By[ind] = Bo[1]*(rf2)-Bo[1]
+        Bz[ind] = Bo[2]*(rf2)-Bo[2]
+
+    r = Utils.mkvc(np.sqrt((x-xc)**2+(y-yc)**2+(z-zc)**2 ))
+    V = 4*np.pi*R**3/3
+    mom = Bot/mu_0*chi/(1+chi/3)*V
+    const = mu_0/(4*np.pi)*mom
+    mdotr = (mx*(x[~ind]-xc)/r[~ind] + my*(y[~ind]-yc)/r[~ind] + mz*(z[~ind]-zc)/r[~ind])
+    Bx[~ind] = const/(r[~ind]**3)*(3*mdotr*(x[~ind]-xc)/r[~ind]-mx)
+    By[~ind] = const/(r[~ind]**3)*(3*mdotr*(y[~ind]-yc)/r[~ind]-my)
+    Bz[~ind] = const/(r[~ind]**3)*(3*mdotr*(z[~ind]-zc)/r[~ind]-mz)
+
+
+    return Bx, By, Bz
+
+
+def IDTtoxyz(Inc, Dec, Btot):
+    """
+        Convert from Inclination, Declination, Total  intensity of earth field to x, y, z
+    """
+    Bx = Btot*np.cos(Inc/180.*np.pi)*np.sin(Dec/180.*np.pi)
+    By = Btot*np.cos(Inc/180.*np.pi)*np.cos(Dec/180.*np.pi)
+    Bz = -Btot*np.sin(Inc/180.*np.pi)
+
+    return np.r_[Bx, By, Bz]
+
+
+def MagSphereFreeSpace(x, y, z, R, xc, yc, zc, chi, Bo):
+    """
+        Computing boundary condition using Congrous sphere method.
+        This is designed for secondary field formulation.
+        >> Input
+        mesh:   Mesh class
+        Bo:     np.array([Box, Boy, Boz]): Primary magnetic flux
+        Chi:    susceptibility at cell volume
+
+        .. math::
+
+            \\vec{B}(r) = \\frac{\mu_0}{4\pi}\\frac{m}{\| \\vec{r}-\\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
+
+    """
+    if (~np.size(x)==np.size(y)==np.size(z)):
+        print "Specify same size of x, y, z"
+        return
+
+    x = Utils.mkvc(x)
+    y = Utils.mkvc(y)
+    z = Utils.mkvc(z)
+
+    nobs = len(x)
+
+    Bot = np.sqrt(sum(Bo**2))
+
+    mx = np.ones([nobs]) * Bo[0,0] *  R**3 / 3. * chi
+    my = np.ones([nobs]) * Bo[0,1] *  R**3 / 3. * chi
+    mz = np.ones([nobs]) * Bo[0,2] *  R**3 / 3. * chi
+
+    M = np.c_[mx, my, mz]
+
+    rx = (x - xc)
+    ry = (y - yc)
+    rz = (zc - z)
+
+    rvec = np.c_[rx, ry, rz]
+    r = np.sqrt((rx)**2+(ry)**2+(rz)**2 )
+
+    B = -Utils.sdiag(1./r**3)*M + Utils.sdiag((3 * np.sum(M*rvec,axis=1))/r**5)*rvec
+
+    Bx = B[:,0]
+    By = B[:,1]
+    Bz = B[:,2]
+
+    return Bx, By, Bz
+
+if __name__ == '__main__':
+
+    hxind = [(0,25,1.3),(21, 12.5),(0,25,1.3)]
+    hyind = [(0,25,1.3),(21, 12.5),(0,25,1.3)]
+    hzind = [(0,25,1.3),(20, 12.5),(0,25,1.3)]
+    # hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)
+    M3 = Mesh.TensorMesh([hxind, hyind, hzind], "CCC")
+    indxd, indxu, indyd, indyu, indzd, indzu = M3.faceBoundaryInd
+    mu0 = 4*np.pi*1e-7
+    chibkg = 0.
+    chiblk = 0.01
+    chi = np.ones(M3.nC)*chibkg
+    sph_ind = spheremodel(M3, 0, 0, 0, 100)
+    chi[sph_ind] = chiblk
+    mu = (1.+chi)*mu0
+    Bbc, const = CongruousMagBC(M3, np.array([1., 0., 0.]), chi)
+
+    flag = 'secondary'
+    Box = 1.
+    H0 = Box/mu_0
+    Bbcxx, Bbcxy, Bbcxz  = MagSphereAnaFun(M3.gridFx[(indxd|indxu),0], M3.gridFx[(indxd|indxu),1], M3.gridFx[(indxd|indxu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
+    Bbcyx, Bbcyy, Bbcyz  = MagSphereAnaFun(M3.gridFy[(indyd|indyu),0], M3.gridFy[(indyd|indyu),1], M3.gridFy[(indyd|indyu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
+    Bbczx, Bbczy, Bbczz  = MagSphereAnaFun(M3.gridFz[(indzd|indzu),0], M3.gridFz[(indzd|indzu),1], M3.gridFz[(indzd|indzu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
+    Bbc_ana = np.r_[Bbcxx, Bbcyy, Bbczz]
+
+    # fig, ax = plt.subplots(1,1, figsize = (10, 10))
+    # ax.plot(Bbc_ana)
+    # ax.plot(Bbc)
+    # plt.show()
+    err = np.linalg.norm(Bbc-Bbc_ana)/np.linalg.norm(Bbc_ana)
+
+    if err < 0.1:
+        print 'Mag Boundary computation is valid, err = ', err
+    else:
+        print 'Mag Boundary computation is wrong!!, err = ', err
+    pass

SimPEG/PF/MagneticsDriver.py

@@ -0,0 +1,334 @@
+import re, os
+from SimPEG import Mesh, np, Utils
+import BaseMag, Magnetics
+
+class MagneticsDriver_Inv(object):
+    """docstring for MagneticsDriver_Inv"""
+
+    def __init__(self, input_file=None):
+        if input_file is not None:
+            self.basePath = os.path.sep.join(input_file.split(os.path.sep)[:-1])
+            if len(self.basePath) > 0:
+                self.basePath += os.path.sep
+            self.readDriverFile(input_file.split(os.path.sep)[-1])
+
+
+    def readDriverFile(self, input_file):
+        """
+            Read input files for forward modeling MAG data with integral form
+            INPUT:
+            input_file: File name containing the forward parameter
+
+            OUTPUT:
+            mshfile
+            obsfile
+            topofile
+            start model
+            ref model
+            mag model
+            weightfile
+            chi_target
+            as, ax ,ay, az
+            upper, lower bounds
+            lp, lqx, lqy, lqz
+
+            # All files should be in the working directory, otherwise the path must
+            # be specified.
+
+        """
+
+
+        fid = open(self.basePath + input_file,'r')
+
+        # Line 1
+        line = fid.readline()
+        l_input  = line.split('!')
+        mshfile = l_input[0].rstrip()
+
+        # Line 2
+        line = fid.readline()
+        l_input  = line.split('!')
+        obsfile = l_input[0].rstrip()
+
+        # Line 3
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input=='null':
+            topofile = []
+
+        else:
+            topofile = l_input[0].rstrip()
+
+
+        # Line 4
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            mstart = float(l_input[1])
+
+        else:
+            mstart = l_input[0].rstrip()
+
+        # Line 5
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            mref = float(l_input[1])
+
+        else:
+            mref = l_input[0].rstrip()
+
+        # Line 6
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            staticInput = float(l_input[1])
+
+        elif l_input[0]=='DEFAULT':
+            staticInput = None
+
+        else:
+            staticInput = l_input[0].rstrip()
+
+
+        # Line 7
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input=='DEFAULT':
+            magfile = []
+
+        else:
+            magfile = l_input[0].rstrip()
+
+        # Line 8
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input=='DEFAULT':
+            wgtfile = []
+
+        else:
+            wgtfile = l_input[0].rstrip()
+
+        # Line 9
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        chi = float(l_input[0])
+
+        # Line 10
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        val = np.array(l_input[0:4])
+        alphas = val.astype(np.float)
+
+        # Line 11
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:3])
+            bounds = val.astype(np.float)
+
+        else:
+            bounds = l_input[0].rstrip()
+
+        # Line 12
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:6])
+            lpnorms = val.astype(np.float)
+
+        else:
+            lpnorms = l_input[0].rstrip()
+
+        # Line 13
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:3])
+            eps = val.astype(np.float)
+
+        else:
+            eps = [None,None]
+
+        self.mshfile  = mshfile
+        self.obsfile  = obsfile
+        self.topofile = topofile
+        self.mstart   = mstart
+        self._mrefInput = mref
+        self._staticInput = staticInput
+        self.magfile  = magfile
+        self.wgtfile  = wgtfile
+        self.chi      = chi
+        self.alphas   = alphas
+        self.bounds   = bounds
+        self.lpnorms  = lpnorms
+        self.eps      = eps
+
+    @property
+    def mesh(self):
+        if getattr(self, '_mesh', None) is None:
+            self._mesh = Mesh.TensorMesh.readUBC(self.basePath + self.mshfile)
+        return self._mesh
+
+    @property
+    def survey(self):
+        if getattr(self, '_survey', None) is None:
+            self._survey = self.readMagneticsObservations(self.obsfile)
+        return self._survey
+
+    @property
+    def activeCells(self):
+        if getattr(self, '_activeCells', None) is None:
+            if self.topofile == 'null':
+                self._activeCells = np.arange(self.mesh.nC)
+            else:
+                topo = np.genfromtxt(self.basePath + self.topofile, skip_header=1)
+                # Find the active cells
+                active = Utils.surface2ind_topo(self.mesh,topo,'N')
+                inds = np.asarray([inds for inds, elem in enumerate(active, 1) if elem], dtype = int) - 1
+                self._activeCells = inds
+
+        return self._activeCells
+
+    @property
+    def staticCells(self):
+        if getattr(self, '_staticCells', None) is None:
+
+            if getattr(self, '_staticInput', None) is None:
+                # All cells are dynamic: 1's
+                self._dynamicCells = np.arange(len(self.m0))
+                self._staticCells = []
+
+            # Cells with specific value are static: 0's
+            else:
+                if isinstance(self._staticInput, float):
+                    staticCells = self.m0 == self._staticInput
+
+                else:
+                    # Read from file active cells with 0:air, 1:dynamic, -1 static
+                    staticCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
+                    staticCells = staticCells[self.activeCells] == -1
+
+                inds = np.asarray([inds for inds, elem in enumerate(staticCells, 1) if elem], dtype = int) - 1
+                self._staticCells = inds
+
+        return self._staticCells
+
+    @property
+    def dynamicCells(self):
+        if getattr(self, '_dynamicCells', None) is None:
+
+            if getattr(self, '_staticInput', None) is None:
+                # All cells are dynamic: 1's
+                self._dynamicCells = np.arange(len(self.m0))
+
+            # Cells with specific value are static: 0's
+            else:
+                if isinstance(self._staticInput, float):
+                    dynamicCells = self.m0 != self._staticInput
+
+                else:
+                    # Read from file active cells with 0:air, 1:dynamic, -1 static
+                    dynamicCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
+                    dynamicCells = dynamicCells[self.activeCells] == 1
+
+                inds = np.asarray([inds for inds, elem in enumerate(dynamicCells, 1) if elem], dtype = int) - 1
+                self._dynamicCells = inds
+
+        return self._dynamicCells
+
+    @property
+    def nC(self):
+        if getattr(self, '_nC', None) is None:
+            self._nC = len(self.activeCells)
+        return self._nC
+
+    @property
+    def m0(self):
+        if getattr(self, '_m0', None) is None:
+            if isinstance(self.mstart, float):
+                self._m0 = np.ones(self.nC) * self.mstart
+            else:
+                self._m0 = Mesh.TensorMesh.readModelUBC(self.mesh,self.basePath + self.mstart)
+                self._m0 = self._m0[self.activeCells]
+
+        return self._m0
+
+    @property
+    def mref(self):
+        if getattr(self, '_mref', None) is None:
+            if isinstance(self._mrefInput, float):
+                self._mref = np.ones(self.nC) * self._mrefInput
+            else:
+                self._mref = Mesh.TensorMesh.readModelUBC(self.mesh,self.basePath + self._mrefInput)
+                self._mref = self._mref[self.activeCells]
+        return self._mref
+
+
+    @property
+    def magnetizationModel(self):
+        """
+            magnetization vector
+        """
+
+        if self.magfile == 'DEFAULT':
+            return Magnetics.dipazm_2_xyz(np.ones(self.nC) * self.survey.srcField.param[1], np.ones(self.nC) * self.survey.srcField.param[2])
+
+        else:
+            raise NotImplementedError("this will require you to read in a three column vector model")
+            self._mref = Utils.meshutils.readUBCTensorModel(self.basePath + self._mrefInput, self.mesh)
+            return np.genfromtxt(self.magfile,delimiter=' \n',dtype=np.str,comments='!')
+
+    def readMagneticsObservations(self, obs_file):
+        """
+            Read and write UBC mag file format
+
+            INPUT:
+            :param fileName, path to the UBC obs mag file
+
+            OUTPUT:
+            :param survey
+            :param M, magnetization orentiaton (MI, MD)
+        """
+
+        fid = open(self.basePath + obs_file,'r')
+
+        # First line has the inclination,declination and amplitude of B0
+        line = fid.readline()
+        B = np.array(line.split(),dtype=float)
+
+        # Second line has the magnetization orientation and a flag
+        line = fid.readline()
+        M = np.array(line.split(),dtype=float)
+
+        # Third line has the number of rows
+        line = fid.readline()
+        ndat = np.array(line.split(),dtype=int)
+
+        # Pre-allocate space for obsx, obsy, obsz, data, uncert
+        line = fid.readline()
+        temp = np.array(line.split(),dtype=float)
+
+        d  = np.zeros(ndat, dtype=float)
+        wd = np.zeros(ndat, dtype=float)
+        locXYZ = np.zeros( (ndat,3), dtype=float)
+
+        for ii in range(ndat):
+
+            temp = np.array(line.split(),dtype=float)
+            locXYZ[ii,:] = temp[:3]
+
+            if len(temp) > 3:
+                d[ii] = temp[3]
+
+                if len(temp)==5:
+                    wd[ii] = temp[4]
+
+            line = fid.readline()
+
+        rxLoc = BaseMag.RxObs(locXYZ)
+        srcField = BaseMag.SrcField([rxLoc],param=(B[2],B[0],B[1]))
+        survey = BaseMag.LinearSurvey(srcField)
+        survey.dobs =  d
+        survey.std =  wd
+        return survey

SimPEG/PF/__init__.py

@@ -0,0 +1,7 @@
+import MagAnalytics
+import BaseMag
+import Magnetics
+import BaseGrav
+import Gravity
+import MagneticsDriver
+import GravityDriver