From 3b97872601fee861a7974ed4dcf3604fb0f50a97 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Sun, 23 Feb 2014 13:18:02 -0800
Subject: [PATCH] Initial work on bringing mag into simpeg framework

---
 simpegPF/BaseMag.py      | 36 +++++++++++++++++++++
 simpegPF/MagAnalytics.py |  2 +-
 simpegPF/Magnetics.py    | 70 ++++++++++++++++++++++++++++++++++++----
 simpegPF/__init__.py     |  4 ++-
 4 files changed, 103 insertions(+), 9 deletions(-)
 create mode 100644 simpegPF/BaseMag.py

diff --git a/simpegPF/BaseMag.py b/simpegPF/BaseMag.py
new file mode 100644
index 00000000..fdd4aea1
--- /dev/null
+++ b/simpegPF/BaseMag.py
@@ -0,0 +1,36 @@
+from SimPEG import Model, Data, np, sp
+from scipy.constants import mu_0
+
+
+class BaseMagData(Data.BaseData):
+    """Base Magnetics Data"""
+
+    rxLoc = None #: receiver locations
+    rxType = None #: receiver type
+
+    def __init__(self, **kwargs):
+        Data.BaseData.__init__(**kwargs)
+
+
+    #TODO: change to inc, dec, intensity
+    def setBackgroundField(self, x, y, z):
+        # Primary field in x-direction (background)
+        self.B0 = np.r_[1.,0,0]
+
+
+class BaseMagModel(Model.BaseModel):
+    """BaseMagModel"""
+
+    def __init__(self, mesh, **kwargs):
+        Model.BaseModel.__init__(mesh, **kwargs)
+
+    def transform(self, m, asMu=True):
+        if asMu:
+            return mu_0*(1 + m)
+        return m
+
+    def transformDeriv(self, m, asMu=True):
+        if asMu:
+            return mu_0*sp.identity(self.nP)
+        return sp.identity(self.nP)
+
diff --git a/simpegPF/MagAnalytics.py b/simpegPF/MagAnalytics.py
index 0504ed45..08fad64e 100644
--- a/simpegPF/MagAnalytics.py
+++ b/simpegPF/MagAnalytics.py
@@ -87,7 +87,7 @@ def CongruousMagBC(mesh, Bo, chi):
     ind = chi > 0.
     V = mesh.vol[ind].sum()
 
-    gamma = 1/V*(chi*mesh.vol).sum()
+    gamma = 1/V*(chi*mesh.vol).sum() # like a mass!
 
     Bot = np.sqrt(sum(Bo**2))
     mx = Bo[0]/Bot
diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py
index 25d375f3..56106db9 100644
--- a/simpegPF/Magnetics.py
+++ b/simpegPF/Magnetics.py
@@ -1,11 +1,62 @@
-from SimPEG import *
+from SimPEG import Problem, Utils
+import BaseMag
+from scipy.constants import mu_0
 
 
-class Magnetics(object):
-    """docstring for Magnetics"""
-    def __init__(self, arg):
-        super(Magnetics, self).__init__()
-        self.arg = arg
+
+class MagneticsDiffSecondary(Problem.BaseProblem):
+    """Secondary field approach using differential equations!"""
+
+    dataPair = BaseMag.BaseMagData
+    modelPair = BaseMag.BaseMagModel
+
+    def __init__(self, mesh, model, **kwargs):
+        Problem.BaseProblem.__init__(mesh, model, **kwargs)
+
+        self._Pbc, self._Pin, self._Pout = /
+            self.mesh.getBCProjWF('neumann', discretization='CC')
+
+        Dface = self.mesh.faceDiv
+        Mc = Utils.sdiag(self.mesh.vol)
+        self._Div = Mc*Dface*self.Pin.T*self.Pin
+
+    @property
+    def Pbc(self): return self._Pbc
+
+    @property
+    def Pin(self): return self._Pin
+
+    @property
+    def Pout(self): return self._Pout
+
+    @property
+    def Div(self): return self._Div
+
+    @property
+    def MfMuI(self): return self._MfMuI
+
+    @property
+    def Mfmu0(self): return self._Mfmu0
+
+    def makeMassMatrices(self, m):
+        mu = self.model.transform(m)
+        Mfmui = self.mesh.getFaceInnerProduct(1./mu)
+        #TODO: this will break if tensor mu
+        self._MfmuI = Utils.sdiag(1./Mfmui.diagonal())
+        self._Mfmu0 = self.mesh.getFaceInnerProduct(1/mu_0)
+
+    def getRHS(self):
+        b0 = self.data.B0
+        B0 = np.r_[b0[0]*np.ones(M3.nFx),
+                   b0[1]*np.ones(M3.nFy),
+                   b0[2]*np.ones(M3.nFz)]
+
+        Dface = self.mesh.faceDiv
+        Mc = Utils.sdiag(self.mesh.vol)
+
+        Bbc = CongruousMagBC(M3, np.array([Box, Boy, Boz]), chi)
+
+        rhs = -self.Div*self.MfmuI*self.Mfmu0*B0 + self.Div*B0 - Mc*Dface*self.Pout.T*Bbc
 
     def getA(self, m):
         """
@@ -20,7 +71,12 @@ class Magnetics(object):
 
 
         """
-        return self.mesh.faceDiv
+        return -self.Div*self.MfmuI*self.Div.T
 
 
+    def fields(self, m):
+        self.makeMassMatrices(m)
+        # F = self.getInitialFields()
+        # return self.forward(m, self.getRHS, self.calcFields, F=F)
+
 
diff --git a/simpegPF/__init__.py b/simpegPF/__init__.py
index 0bc4e552..03bd44cc 100644
--- a/simpegPF/__init__.py
+++ b/simpegPF/__init__.py
@@ -1 +1,3 @@
-# Init!
+import MagAnalytics
+import MagData
+import Magnetics