diff --git a/simpegEM/Base.py b/simpegEM/Base.py index d9ea7984..3bf0dec7 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -39,101 +39,18 @@ class BaseEMProblem(Problem.BaseProblem): self.__makeASymmetric = True return self.__makeASymmetric - #################################################### - # Phys Props - #################################################### - - # Mu - # @property - # def mu(self): - # if getattr(self, '_mu', None) is None: - # # if getattr(self, '_mui', None) is not None: - # # self._mu = sel - # self._mu = mu_0 - # return self._mu - # @mu.setter - # def mu(self, value): - # if getattr(self, '_MfMui', None) is not None: - # del self._MfMui - # if getattr(self, '_MfMuiI', None) is not None: - # del self._MfMuiI - # if getattr(self, '_MeMu', None) is not None: - # del delf._MeMu - # if getattr(self, '_MeMuI', None) is not None: - # del self._MeMuI - # self._mu = value - - # TODO: hardcoded to assume diagonal mu - # @property - # def mui(self): - # if getattr(self, '_mui', None) is None: - # self._mui = 1./mu_0 - # return self._mui - # @mui.setter - # def mui(self, value): - # if getattr(self, '_MfMui', None) is not None: - # del self._MfMui - # if getattr(self, '_MfMuiI', None) is not None: - # del self._MfMuiI - # if getattr(self, '_MeMu', None) is not None: - # del delf._MeMu - # if getattr(self, '_MeMuI', None) is not None: - # del self._MeMuI - # self._mui = value - - # Sigma - # @property - # def sigma(self): - # if getattr(self, '_sigma', None) is None: - # self._sigma = self.curModel.transform - # return self._sigma - # @sigma.setter - # def sigma(self, value): - # if getattr(self, '_MeSigma', None) is not None: - # del self._MeSigma - # if getattr(self, '_MeSigmaI', None) is not None: - # del self._MeSigmaI - # if getattr(self, '_MfSigmai', None) is not None: - # del delf._MfSigmai - # if getattr(self, '_MfSigmaiI', None) is not None: - # del self._MfSigmaiI - # self._sigma = value - - # def dsigma_dm(self): - # return self.curModel.transformDeriv - - - # TODO: hardcoded to assume diagonal sigma - # @property - # def sigmai(self): - # if getattr(self, '_sigmai', None) is None: - # self._sigmai = 1./self.curModel.transform - # return self._sigmai - # @sigmai.setter - # def sigmai(self, value): - # if getattr(self, '_MeSigma', None) is not None: - # del self._MeSigma - # if getattr(self, '_MeSigmaI', None) is not None: - # del self._MeSigmaI - # if getattr(self, '_MfSigmai', None) is not None: - # del delf._MfSigmai - # if getattr(self, '_MfSigmaiI', None) is not None: - # del self._MfSigmaiI - # self._sigma = value - #################################################### # Mass Matrices #################################################### - # TODO: Link to EMPropMap - # if Prop - # deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai','_MfSigmaiI'] @property def deleteTheseOnModelUpdate(self): toDelete = [] - if self.mapping.sigmaMap is not None: - toDelete += ['_MeSigma', '_MeSigmaI','_MfSigmai','_MfSigmaiI'] + if self.mapping.sigmaMap is not None or self.mapping.rhoMap is not None: + toDelete += ['_MeSigma', '_MeSigmaI','_MfRho','_MfRhoI'] + if self.mapping.muMap is not None or self.mapping.muiMap is not None: + toDelete += ['_MeMu', '_MeMuI','_MfMui','_MfMuiI'] return toDelete @property @@ -194,20 +111,20 @@ class BaseEMProblem(Problem.BaseProblem): # def dMeSigmaI_dsigma(self,u) @property - def MfSigmai(self): - if getattr(self, '_MfSigmai', None) is None: - self._MfSigmai = self.mesh.getFaceInnerProduct(self.curModel.rho) - return self._MfSigmai + def MfRho(self): + if getattr(self, '_MfRho', None) is None: + self._MfRho = self.mesh.getFaceInnerProduct(self.curModel.rho) + return self._MfRho - # def dMfSigmai_dsigmai(self,u) + # def dMfRho_dsigmai(self,u) @property - def MfSigmaiI(self): - if getattr(self, '_MfSigmaiI', None) is None: - self._MfSigmaiI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True) - return self._MfSigmaiI + def MfRhoI(self): + if getattr(self, '_MfRhoI', None) is None: + self._MfRhoI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True) + return self._MfRhoI - # def dMfSigmaiI(self,u) + # def dMfRhoI(self,u) #################################################### diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index e83fe0e0..7ba70da2 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -327,14 +327,14 @@ class ProblemFDEM_j(BaseFDEMProblem): """ MeMuI = self.MeMuI - MfSigmai = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - A = C * MeMuI * C.T * MfSigmai + iomega + A = C * MeMuI * C.T * MfRho + iomega if self._makeASymmetric is True: - return MfSigmai.T*A + return MfRho.T*A return A @@ -347,7 +347,7 @@ class ProblemFDEM_j(BaseFDEMProblem): """ MeMuI = self.MeMuI - MfSigmai = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl sigi = self.sigmai dsig_dm = self.curModel.transformDeriv @@ -356,11 +356,11 @@ class ProblemFDEM_j(BaseFDEMProblem): if adjoint: if self._makeASymmetric is True: - v = MfSigmai * v + v = MfRho * v return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) ) if self._makeASymmetric is True: - return MfSigmai.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) + return MfRho.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) @@ -378,8 +378,8 @@ class ProblemFDEM_j(BaseFDEMProblem): RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e if self._makeASymmetric is True: - MfSigmai = self.MfSigmai - return MfSigmai.T*RHS + MfRho = self.MfRho + return MfRho.T*RHS return RHS @@ -429,10 +429,10 @@ class ProblemFDEM_h(BaseFDEMProblem): """ MeMu = self.MeMu - MfSigmai = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl - return C.T * MfSigmai * C + 1j*omega(freq)*MeMu + return C.T * MfRho * C + 1j*omega(freq)*MeMu def getADeriv(self, freq, u, v, adjoint=False): @@ -458,9 +458,9 @@ class ProblemFDEM_h(BaseFDEMProblem): S_m, S_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl - MfSigmai = self.MfSigmai + MfRho = self.MfRho - RHS = S_m + C.T * ( MfSigmai * S_e ) + RHS = S_m + C.T * ( MfRho * S_e ) return RHS diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index e12752b6..9e93d3d5 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -114,7 +114,7 @@ class FieldsFDEM_j(FieldsFDEM): def startup(self): self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeMuI = self.survey.prob.MeMuI - self._MfSigmai = self.survey.prob.MfSigmai + self._MfRho = self.survey.prob.MfRho self._curModel = self.survey.prob.curModel def _j(self, j_sol, srcList): @@ -128,9 +128,9 @@ class FieldsFDEM_j(FieldsFDEM): def _h(self, j_sol, srcList): MeMuI = self._MeMuI C = self._edgeCurl - MfSigmai = self._MfSigmai + MfRho = self._MfRho - h = MeMuI * (C.T * (MfSigmai * j_sol) ) + h = MeMuI * (C.T * (MfRho * j_sol) ) for i, src in enumerate(srcList): h[:,i] *= -1./(1j*omega(src.freq)) @@ -151,7 +151,7 @@ class FieldsFDEM_j(FieldsFDEM): dsig_dm = self._curModel.transformDeriv dsigi_dsig = -Utils.sdiag(sigi)**2 dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) - sigi = self._MfSigmai + sigi = self._MfRho S_mDeriv,_ = src.getSourceDeriv(self.survey.prob, v, adjoint) @@ -177,7 +177,7 @@ class FieldsFDEM_h(FieldsFDEM): def startup(self): self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeMuI = self.survey.prob.MeMuI - self._MfSigmai = self.survey.prob.MfSigmai + self._MfRho = self.survey.prob.MfRho def _h(self, h_sol, srcList): h = h_sol @@ -215,11 +215,11 @@ class FieldsFDEM_h(FieldsFDEM): # elif fieldType == 'h': # MeMuI = self._MeMuI # C = self.mesh.edgeCurl - # MfSigmai = self._MfSigmai + # MfRho = self._MfRho # if not adjoint: - # h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfSigmai * j ) ) + # h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfRho * j ) ) # else: - # h = -(1./(1j*omega(freq))) * MfSigmai.T * ( C * ( MeMuI.T * j ) ) + # h = -(1./(1j*omega(freq))) * MfRho.T * ( C * ( MeMuI.T * j ) ) # return h # raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) @@ -235,7 +235,7 @@ class FieldsFDEM_h(FieldsFDEM): # dsig_dm = self._curModel.transformDeriv # dsigi_dsig = -Utils.sdiag(sigi)**2 # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) - # sigi = self._MfSigmai + # sigi = self._MfRho # if not adjoint: # return -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) # else: diff --git a/simpegEM/Utils/EMUtils.py b/simpegEM/Utils/EMUtils.py index 5cdaf150..916d6ae2 100644 --- a/simpegEM/Utils/EMUtils.py +++ b/simpegEM/Utils/EMUtils.py @@ -16,7 +16,7 @@ def e_from_j(prob,j): if eqLocs is 'FE': MSigmaI = prob.MeSigmaI elif eqLocs is 'EF': - MSigmaI = prob.MfSigmai + MSigmaI = prob.MfRho return MSigmaI*j def j_from_e(prob,e): @@ -24,13 +24,13 @@ def j_from_e(prob,e): if eqLocs is 'FE': MSigma = prob.MeSigma elif eqLocs is 'EF': - MSigma = prob.MfSigma + MSigma = prob.MfRhoI return MSigma*e def b_from_h(prob,h): eqLocs = prob._eqLocs if eqLocs is 'FE': - MMu = prob.MfMu + MMu = prob.MfMuiI elif eqLocs is 'EF': MMu = prob.MeMu return MMu*h