1. Add distributed source for nodal discretization

2. Add Analytic tests
3. Fix simple bug in PlotSlice for nodal variable
4. Add more analytic function (sphere)

SimPEG/EM/Analytics/DC.py

@@ -1,7 +1,7 @@
 import numpy as np
 from scipy.constants import mu_0, pi
 
-def DCAnalytic(txloc, rxlocs, sigma, flag="wholespace"):
+def DCAnalyticHalf(txloc, rxlocs, sigma, flag="wholespace"):
     """
         Analytic solution for electric potential from a postive pole
 
@@ -32,3 +32,82 @@ def DCAnalytic(txloc, rxlocs, sigma, flag="wholespace"):
 
     return phi
 
+deg2rad  = lambda deg: deg/180.*np.pi
+rad2deg  = lambda rad: rad*180./np.pi
+
+def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
+                 flag = "sec", order=12):
+# def DCSpherePointCurrent(txloc, rxloc, xc, radius, rho, rho1, \
+#                  flag = "sec", order=12):
+    """
+
+        Parameters:
+
+            txloc (array) : current electrode location (x,y,z)
+            xc (float)    : x center of depressed sphere
+            rxloc (array) : electrode locations
+                              (Nx3 array, # of electrodes)
+            radius (float): radius of the sphere (m)
+            rho (float)   : resistivity of the background (ohm-m)
+            rho1 (float)  : resistivity of the sphere
+            flag (string) : "sec", "total", "prim"
+                              (default="sec")
+                              "sec": secondary potential only due to sphere
+                              "prim": primary potential from the point source
+                              "total": "sec"+"prim"
+            order (float) : 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 flag =="prim":
+        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 flag == "sec":
+        return out-prim
+    elif flag == "total":
+        return 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,4 +1,4 @@
 from TDEM import hzAnalyticDipoleT
 from FDEM import hzAnalyticDipoleF
 from FDEMcasing import *
-from DC import DCAnalytic
+from DC import DCAnalyticHalf, DCAnalyticSphere

SimPEG/EM/Static/DC/SrcDC.py

@@ -1,6 +1,6 @@
 import SimPEG
 # from SimPEG.EM.Base import BaseEMSurvey
-from SimPEG.Utils import Zero, closestPoints
+from SimPEG.Utils import Zero, closestPoints, mkvc
 import numpy as np
 
 class BaseSrc(SimPEG.Survey.BaseSrc):
@@ -27,13 +27,13 @@ class Dipole(BaseSrc):
 
     def eval(self, prob):
         if prob._formulation == 'HJ':
-            inds = closestPoints(prob.mesh, self.loc)
+            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':
-            inds = closestPoints(prob.mesh, self.loc)
-            q = np.zeros(prob.mesh.nN)
-            q[inds] = self.current * np.r_[1., -1.]
+            qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense()
+            qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense()
+            q = mkvc(qa+qb)
         return q
 
     # def bc_contribution

SimPEG/Mesh/View.py

@@ -206,7 +206,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)

tests/em/static/test_DC_analytic.py

@@ -0,0 +1,68 @@
+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, flag="halfspace")
+        phiB = EM.Analytics.DCAnalyticHalf(Bloc, [M,N], 1e-2, flag="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_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
+        else:
+            passed = False
+        self.assertTrue(passed)
+
+    def test_CC(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_CC is pased"
+        else:
+            passed = False
+        self.assertTrue(passed)
+
+if __name__ == '__main__':
+    unittest.main()
+