From ed1457df86a4d22e30c312bb48525686ce1e7452 Mon Sep 17 00:00:00 2001
From: Lindsey <lindseyheagy@gmail.com>
Date: Mon, 23 Feb 2015 17:42:25 -0800
Subject: [PATCH] start of j implementation, UNTESTED

---
 simpegEM/Base.py          | 17 +++++++++-
 simpegEM/FDEM/FDEM.py     | 69 +++++++++++++++++++++++++++++++++++++++
 simpegEM/FDEM/__init__.py |  2 +-
 3 files changed, 86 insertions(+), 2 deletions(-)

diff --git a/simpegEM/Base.py b/simpegEM/Base.py
index 373cee75..da7e4a86 100644
--- a/simpegEM/Base.py
+++ b/simpegEM/Base.py
@@ -28,6 +28,13 @@ class BaseEMProblem(Problem.BaseProblem):
             self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0)
         return self._MfMui
 
+    @property
+    def MeMuI(self):
+        #TODO: assuming constant mu
+        if getattr(self, '_MeMuI', None) is None:
+            self._MeMuI = self.mesh.getEdgeInnerProduct(1/mu_0)
+        return self._MeMuI
+
     @property
     def Me(self):
         if getattr(self, '_Me', None) is None:
@@ -50,7 +57,15 @@ class BaseEMProblem(Problem.BaseProblem):
             self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True)
         return self._MeSigmaI
 
-    deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI']
+    @property
+    def MfSigmaI(self):
+        #TODO: hardcoded to sigma as the model
+        if getattr(self, '_MfSigmaI', None) is None:
+            sigma = self.curModel.transform
+            self._MfSigmaI = self.mesh.getFaceInnerProduct(sigma, invMat=True)
+        return self._MfSigmaI
+
+    deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmaI']
 
     def fields(self, m):
         self.curModel = m
diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py
index c4137486..3f0af1fb 100644
--- a/simpegEM/FDEM/FDEM.py
+++ b/simpegEM/FDEM/FDEM.py
@@ -335,3 +335,72 @@ class ProblemFDEM_b(BaseFDEMProblem):
             return None
         raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
 
+
+class ProblemFDEM_j(BaseFDEMProblem):
+    """
+        Using the H-J formulation of Maxwell's equations
+
+        .. math::
+            \\nabla \\times \\sigma^{-1} \\vec{J} + i\\omega\\mu\\vec{H} = 0
+            \\nabla \\times \\vec{H} - \\vec{J} = \\vec{J_s}
+
+        Since \(\\vec{J}\) is a flux and \(\\vec{H}\) is a field, we discretize \(\\vec{J}\) on faces and \(\\vec{H}\) on edges.
+
+        For this implementation, we solve for J, using \( \\vec{H} = - (i\\omega\\mu)^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} \) :
+
+        .. math::
+            \\nabla \\times ( \\mu^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} ) + i\\omega \\vec{J} = - i\\omega\\vec{J_s}
+    """
+
+    solType = 'j'
+
+    def __init__(self, model, **kwargs):
+        BaseFDEMProblem.__init__(self, model, **kwargs)
+
+    def getA(self, freq):
+        """
+            :param float freq: Frequency
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+
+        muI = self.MeMuI
+        sigI = self.MfSigmaI
+        C = self.mesh.edgeCurl
+
+        return C * muI * C.T * sigI + 1j * omega(freq)
+
+    def getADeriv(self, freq, u, v, adjoint=False):
+
+        muI = self.MeMuI
+        C = self.mesh.edgeCurl
+        sig = self.curModel.transform
+        dsig_dm = self.curModel.transformDeriv
+        #TODO: This only works if diagonal (no tensors)...
+        dMfSigmaI_dI = - self.MfSigmaI**2
+
+        dMf_dsig = self.mesh.getFaceInnerProductDeriv(sig)(u)
+
+        if adjoint:
+            raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+
+        return C * ( muI * ( C.T * ( dMfSigmaI_dI * ( dMf_dsig * ( dsig_dm * v ) ) ) ) )
+
+
+        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+
+        
+
+    def getRHS(self, freq):
+        """
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, nTx)
+            :return: RHS
+        """
+        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+
+    def calcFields(self, sol, freq, fieldType, adjoint=False):
+        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+
+    def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
+        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
diff --git a/simpegEM/FDEM/__init__.py b/simpegEM/FDEM/__init__.py
index 70124686..4f7fcf4c 100644
--- a/simpegEM/FDEM/__init__.py
+++ b/simpegEM/FDEM/__init__.py
@@ -1,2 +1,2 @@
 from SurveyFDEM import *
-from FDEM import ProblemFDEM_e, ProblemFDEM_b
+from FDEM import ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j