Extend BaseScr to that Pole and Dipole DC are simply minor extensions…

SimPEG/electromagnetics/static/resistivity/receivers.py

@@ -144,8 +144,8 @@ def __init__(self, locations_m=None, locations_n=None, locations=None, **kwargs)
             )
 
         # if locations_m set, then use locations_m, locations_n
-        if locations_m is not None:
-            if locations_n is None:
+        if locations_m is not None or locations_n is not None:
+            if locations_n is None or locations_m is None:
                 raise ValueError(
                     "For a dipole source both locations_m and locations_n "
                     "must be set"

SimPEG/electromagnetics/static/resistivity/simulation_2d.py

@@ -2,6 +2,8 @@
 from scipy.optimize import minimize
 import warnings
 import properties
+
+import discretize
 from ....utils.code_utils import deprecate_class
 
 from ....utils import mkvc, sdiag, Zero
@@ -60,9 +62,13 @@ def g(k):
 
                 return phi, g
 
-            # find the minimum cell spacing, and the maximum side of the mesh
-            min_r = min(*[np.min(h) for h in self.mesh.h])
-            max_r = max(*[np.sum(h) for h in self.mesh.h])
+            if isinstance(self.mesh, discretize.CurvilinearMesh):
+                min_r = min(self.mesh.edge_lengths)
+                max_r = max(self.mesh.edge_lengths)
+            else:
+                # find the minimum cell spacing, and the maximum side of the mesh
+                min_r = min(*[np.min(h) for h in self.mesh.h])
+                max_r = max(*[np.sum(h) for h in self.mesh.h])
             # generate test points log spaced between these two end members
             rs = np.logspace(np.log10(min_r / 4), np.log10(max_r * 4), 100)
 
@@ -113,7 +119,7 @@ def g(k):
         if miniaturize:
             self._dipoles, self._invs, self._mini_survey = _mini_pole_pole(self.survey)
 
-    def fields(self, m):
+    def fields(self, m=None):
         if self.verbose:
             print(">> Compute fields")
         if m is not None:
@@ -163,7 +169,7 @@ def dpred(self, m=None, f=None):
         for src in survey.source_list:
             for rx in src.receiver_list:
                 d = rx.eval(src, self.mesh, f).dot(weights)
-                temp[count : count + len(d)] = d
+                temp[count: count + len(d)] = d
                 count += len(d)
 
         return self._mini_survey_data(temp)
@@ -223,7 +229,7 @@ def Jvec(self, m, v, f=None):
                     df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False)
                     Jv1_temp = rx.evalDeriv(src, self.mesh, f, df_dm_v)
                     # Trapezoidal intergration
-                    Jv[count : count + len(Jv1_temp)] += weights[iky] * Jv1_temp
+                    Jv[count: count + len(Jv1_temp)] += weights[iky] * Jv1_temp
                     count += len(Jv1_temp)
 
         return self._mini_survey_data(Jv)
@@ -271,7 +277,7 @@ def _Jtvec(self, m, v=None, f=None):
                     df_duT_sum = 0
                     df_dmT_sum = 0
                     for rx in src.receiver_list:
-                        my_v = v[count : count + rx.nD]
+                        my_v = v[count: count + rx.nD]
                         count += rx.nD
                         # wrt f, need possibility wrt m
                         PTv = rx.evalDeriv(src, self.mesh, f, my_v, adjoint=True)
@@ -587,7 +593,7 @@ def setBC(self, ky=None):
                 is_t[:, -1] = True
                 is_t = is_t.reshape(-1, order="F")[is_b]
                 not_top = np.zeros(boundary_faces.shape[0], dtype=bool)
-                not_top[-len(is_t) :] = ~is_t
+                not_top[-len(is_t):] = ~is_t
 
             # use the exponentialy scaled modified bessel function of second kind,
             # (the division will cancel out the scaling)
@@ -742,7 +748,7 @@ def setBC(self, ky=None):
                 is_t[:, -1] = True
                 is_t = is_t.reshape(-1, order="F")[is_b]
                 not_top = np.zeros(boundary_faces.shape[0], dtype=bool)
-                not_top[-len(is_t) :] = ~is_t
+                not_top[-len(is_t):] = ~is_t
 
             # use the exponentiall scaled modified bessel function of second kind,
             # (the division will cancel out the scaling)

SimPEG/electromagnetics/static/resistivity/sources.py

@@ -2,31 +2,90 @@
 import properties
 
 from .... import survey
-from ....utils import Zero, closestPoints
-from ....utils.code_utils import deprecate_property
-
-import warnings
+from ....utils import Zero
 
 
 class BaseSrc(survey.BaseSrc):
     """
     Base DC source
     """
 
-    current = properties.Float("amplitude of the source current", default=1.0)
+    current = properties.List(doc="amplitudes of the source currents", default=[1.0])
+    location = properties.List(
+        "location of the source electrodes",
+        survey.SourceLocationArray("location of electrode"),
+    )
 
     _q = None
 
-    def __init__(self, receiver_list, **kwargs):
+    def __init__(self, receiver_list, location, current=None, **kwargs):
         super(BaseSrc, self).__init__(receiver_list, **kwargs)
+        if type(location) is np.ndarray:
+            if location.ndim == 2:
+                location = list(location)
+            else:
+                location = [location]
+        if current is None:
+            current = self.current
+        if type(current) == float or type(current) == int:
+            current = [float(current)]
+        elif type(current) == np.ndarray:
+            if current.ndim == 2:
+                current = list(current)
+            else:
+                location = [location]
+        if len(current) != 1 and len(current) != len(location):
+            raise ValueError(
+                "Current must be constant or equal to the number of specified source locations."
+            )
+        if type(current) != list:
+            current = np.repeat(current, len(location)).tolist()
+
+        self.current = current
+        self.location = location
 
     def eval(self, sim):
-        raise NotImplementedError
+        if self._q is not None:
+            return self._q
+        else:
+            if sim._formulation == "HJ":
+                inds = sim.mesh.closest_points_index(self.location, grid_loc="CC")
+                self._q = np.zeros(sim.mesh.nC)
+                self._q[inds] = self.current
+            elif sim._formulation == "EB":
+                loc = np.row_stack(self.location)
+                cur = np.asarray(self.current)
+                interpolation_matrix = sim.mesh.get_interpolation_matrix(loc, locType='N').toarray()
+                q = np.sum(cur[:, np.newaxis] * interpolation_matrix, axis=0)
+                self._q = q
+            return self._q
 
     def evalDeriv(self, sim):
         return Zero()
 
 
+class Multipole(BaseSrc):
+    """
+    Generic Multipole Source
+    """
+
+    def __init__(self, receiver_list=[], location=None, **kwargs):
+        super(Multipole, self).__init__(receiver_list, location, **kwargs)
+
+    @property
+    def location_a(self):
+        """Locations of the A electrode"""
+        return self.location
+
+    @property
+    def location_b(self):
+        """Location of the B electrode"""
+        return list(np.tile((len(self.location), 1),
+                            np.full_like(self.location[0], np.nan)
+                            )
+                    )
+
+
 class Dipole(BaseSrc):
     """
     Dipole source
@@ -36,9 +95,6 @@ class Dipole(BaseSrc):
         "location of the source electrodes",
         survey.SourceLocationArray("location of electrode"),
     )
-    loc = deprecate_property(
-        location, "loc", new_name="location", removal_version="0.16.0", error=True
-    )
 
     def __init__(
         self,
@@ -62,6 +118,11 @@ def __init__(
                 "The locationB property has been removed. Please set the "
                 "location_b property instead.",
             )
+        if "current" in kwargs.keys():
+            value = kwargs.pop('current')
+            current = [value, -value]
+        else:
+            current = [1.0, -1.0]
 
         # if location_a set, then use location_a, location_b
         if location_a is not None:
@@ -95,8 +156,7 @@ def __init__(
             )
 
         # instantiate
-        super(Dipole, self).__init__(receiver_list, **kwargs)
-        self.location = location
+        super(Dipole, self).__init__(receiver_list, location, current, **kwargs)
 
     def __repr__(self):
         return (
@@ -113,48 +173,17 @@ def location_b(self):
         """Location of the B-electrode"""
         return self.location[1]
 
-    def eval(self, sim):
-        if self._q is not None:
-            return self._q
-        else:
-            if sim._formulation == "HJ":
-                inds = closestPoints(sim.mesh, self.location, gridLoc="CC")
-                self._q = np.zeros(sim.mesh.nC)
-                self._q[inds] = self.current * np.r_[1.0, -1.0]
-            elif sim._formulation == "EB":
-                qa = sim.mesh.getInterpolationMat(
-                    self.location[0], locType="N"
-                ).toarray()
-                qb = -sim.mesh.getInterpolationMat(
-                    self.location[1], locType="N"
-                ).toarray()
-                self._q = self.current * (qa + qb)
-            return self._q
-
 
 class Pole(BaseSrc):
     def __init__(self, receiver_list=[], location=None, **kwargs):
-        super(Pole, self).__init__(receiver_list, location=location, **kwargs)
-
-    def eval(self, sim):
-        if self._q is not None:
-            return self._q
-        else:
-            if sim._formulation == "HJ":
-                inds = closestPoints(sim.mesh, self.location)
-                self._q = np.zeros(sim.mesh.nC)
-                self._q[inds] = self.current * np.r_[1.0]
-            elif sim._formulation == "EB":
-                q = sim.mesh.getInterpolationMat(self.location, locType="N")
-                self._q = self.current * q.toarray()
-            return self._q
+        super(Pole, self).__init__(receiver_list, location, **kwargs)
 
     @property
     def location_a(self):
         """Locations of the A electrode"""
-        return self.location
+        return self.location[0]
 
     @property
     def location_b(self):
         """Location of the B electrode"""
-        return np.nan * np.ones_like(self.location)
+        return np.nan * np.ones_like(self.location[0])

SimPEG/electromagnetics/static/resistivity/survey.py

@@ -4,11 +4,9 @@
 from __future__ import unicode_literals
 
 import numpy as np
-from scipy.interpolate import interp1d, NearestNDInterpolator
 import properties
 from ....utils.code_utils import deprecate_class, deprecate_property
 
-from ....utils import uniqueRows
 from ....survey import BaseSurvey
 from ..utils import drapeTopotoLoc
 from . import receivers as Rx
@@ -158,7 +156,7 @@ def source_locations(self):
         return [np.vstack(src_a), np.vstack(src_b)]
 
     def set_geometric_factor(
-        self, space_type="half-space", data_type=None, survey_type=None,
+            self, space_type="half-space", data_type=None, survey_type=None,
     ):
         if data_type is not None:
             raise TypeError(
@@ -190,10 +188,10 @@ def _set_abmn_locations(self):
                 # Pole Source
                 if isinstance(source, Src.Pole):
                     locations_a.append(
-                        source.location.reshape([1, -1]).repeat(nRx, axis=0)
+                        source.location[0].reshape([1, -1]).repeat(nRx, axis=0)
                     )
                     locations_b.append(
-                        source.location.reshape([1, -1]).repeat(nRx, axis=0)
+                        source.location[0].reshape([1, -1]).repeat(nRx, axis=0)
                     )
                 # Dipole Source
                 elif isinstance(source, Src.Dipole):
@@ -203,6 +201,15 @@ def _set_abmn_locations(self):
                     locations_b.append(
                         source.location[1].reshape([1, -1]).repeat(nRx, axis=0)
                     )
+                elif isinstance(source, Src.Multipole):
+                    location_as_array = np.asarray(source.location)
+                    location_as_array = np.tile(location_as_array, (nRx, 1))
+                    locations_a.append(
+                        location_as_array
+                    )
+                    locations_b.append(
+                        location_as_array
+                    )
                 # Pole RX
                 if isinstance(rx, Rx.Pole):
                     locations_m.append(rx.locations)
@@ -226,7 +233,7 @@ def getABMN_locations(self):
         )
 
     def drape_electrodes_on_topography(
-        self, mesh, actind, option="top", topography=None, force=False
+            self, mesh, actind, option="top", topography=None, force=False
     ):
         """Shift electrode locations to be on [top] of the active cells."""
         if self.survey_geometry == "surface":
@@ -237,9 +244,9 @@ def drape_electrodes_on_topography(
             unique_electrodes, inv = np.unique(
                 np.vstack((loc_a, loc_b, loc_m, loc_n)), return_inverse=True, axis=0
             )
-            inv_a, inv = inv[: len(loc_a)], inv[len(loc_a) :]
-            inv_b, inv = inv[: len(loc_b)], inv[len(loc_b) :]
-            inv_m, inv_n = inv[: len(loc_m)], inv[len(loc_m) :]
+            inv_a, inv = inv[: len(loc_a)], inv[len(loc_a):]
+            inv_b, inv = inv[: len(loc_b)], inv[len(loc_b):]
+            inv_m, inv_n = inv[: len(loc_m)], inv[len(loc_m):]
 
             electrodes_shifted = drapeTopotoLoc(
                 mesh, unique_electrodes, actind=actind, option=option
@@ -252,8 +259,8 @@ def drape_electrodes_on_topography(
             ind = 0
             for src in self.source_list:
                 a_loc, b_loc = a_shifted[ind], b_shifted[ind]
-                if isinstance(src, Src.Pole):
-                    src.location = a_loc
+                if type(src) is Src.Pole or type(src) is Src.BaseSrc:
+                    src.location = [a_loc]
                 else:
                     src.location = [a_loc, b_loc]
                 for rx in src.receiver_list: