Added writeUBC_DCobs to StaticUtils from DCIPUtils

SimPEG/EM/Static/DC/ProblemDC.py

@@ -35,10 +35,11 @@ class BaseDCProblem(BaseEMProblem):
 
         self.curModel = m
 
-        Jv = self.dataPair(self.survey) #same size as the data
-
+        # 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)
@@ -48,8 +49,10 @@ class BaseDCProblem(BaseEMProblem):
             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)
-        return Utils.mkvc(Jv)
+                # 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:
@@ -64,7 +67,6 @@ class BaseDCProblem(BaseEMProblem):
         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:

SimPEG/EM/Static/IP/ProblemIP.py

@@ -45,7 +45,8 @@ class BaseIPProblem(BaseEMProblem):
 
         self.curModel = m
 
-        Jv = self.dataPair(self.survey) #same size as the data
+        # Jv = self.dataPair(self.survey) #same size as the data
+        Jv = []
 
         A = self.getA()
 
@@ -58,13 +59,16 @@ class BaseIPProblem(BaseEMProblem):
             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[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 -Utils.mkvc(Jv)
+            return -np.hstack(Jv)
         # Conductivity (d u / d log rho)
         if self._formulation is 'HJ':
-            return Utils.mkvc(Jv)
+            # return Utils.mkvc(Jv)
+            return np.hstack(Jv)
 
     def Jtvec(self, m, v, f=None):
         if f is None:

SimPEG/EM/Static/Utils/StaticUtils.py

@@ -315,3 +315,107 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
 
     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/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'