Skip to content

Commit

Permalink
the forward for all fields and fluxes can be computed from any formul…
Browse files Browse the repository at this point in the history
…ation (todo: derive)
  • Loading branch information
lheagy committed Feb 5, 2016
1 parent 8a7e72f commit 73edef5
Show file tree
Hide file tree
Showing 4 changed files with 42 additions and 63 deletions.
10 changes: 5 additions & 5 deletions SimPEG/EM/FDEM/FieldsFDEM.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,8 +20,8 @@ class Fields_e(Fields):
'b' : ['eSolution','F','_b'],
'bPrimary' : ['eSolution','F','_bPrimary'],
'bSecondary' : ['eSolution','F','_bSecondary'],
'j' : ['eSolution','CC','_j'],
'h' : ['eSolution','CC','_h'],
'j' : ['eSolution','CCV','_j'],
'h' : ['eSolution','CCV','_h'],
}

def __init__(self,mesh,survey,**kwargs):
Expand Down Expand Up @@ -189,7 +189,7 @@ def _GLoc(self,fieldType):
elif fieldType == 'b':
return 'F'
elif (fieldType == 'h') or (fieldType == 'j'):
return'CC'
return'CCV'
else:
raise Exception('Field type must be e, b, h, j')

Expand Down Expand Up @@ -325,7 +325,7 @@ def _GLoc(self,fieldType):
elif fieldType == 'j':
return 'F'
elif (fieldType == 'e') or (fieldType == 'b'):
return 'CC'
return 'CCV'
else:
raise Exception('Field type must be e, b, h, j')

Expand Down Expand Up @@ -467,7 +467,7 @@ def _GLoc(self,fieldType):
elif fieldType == 'j':
return 'F'
elif (fieldType == 'e') or (fieldType == 'b'):
return 'CC'
return 'CCV'
else:
raise Exception('Field type must be e, b, h, j')

Expand Down
33 changes: 0 additions & 33 deletions SimPEG/EM/FDEM/SurveyFDEM.py
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,6 @@ def projComp(self):

def projGLoc(self, u):
"""Grid Location projection (e.g. Ex Fy ...)"""
print 'here', u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]

def projectFields(self, src, mesh, u):
Expand All @@ -72,30 +71,6 @@ def projectFields(self, src, mesh, u):
# get the real or imag component
real_or_imag = self.projComp
u_part = getattr(u_part_complex, real_or_imag)



# if projGLoc == 'CC':
# P = self.getP(mesh, projGLoc)
# Z = 0.*P
# if mesh.dim == 3:
# if mesh._meshType == 'CYL' and mesh.isSymmetric and u_part.size > mesh.nC: # TODO: there must be a better way to do this!
# if self.knownRxTypes[self.rxType][1] == 'x':
# P = sp.hstack([P,Z])
# elif self.knownRxTypes[self.rxType][1] == 'z':
# P = sp.hstack([Z,P])
# elif self.knownRxTypes[self.rxType][1] == 'y':
# raise Exception('Symmetric CylMesh does not support y interpolation, as this variable does not exist.')
# else:
# if self.knownRxTypes[self.rxType][1] == 'x':
# P = sp.hstack([P,Z,Z])
# elif self.knownRxTypes[self.rxType][1] == 'y':
# P = sp.hstack([Z,P,Z])
# elif self.knownRxTypes[self.rxType][1] == 'z':
# P = sp.hstack([Z,Z,P])
# else:
# projGLoc += self.knownRxTypes[self.rxType][1]
# P = self.getP(mesh, projGLoc)

return P*u_part

Expand All @@ -106,12 +81,6 @@ def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False):
print self.knownRxTypes[self.rxType][:2], 'Deriv', projGLoc
projGLoc += self.knownRxTypes[self.rxType][1]

# if projGLoc = 'CC':
# P = self.getP(mesh)
# if sel

# else projGLoc != 'CC':
# projGLoc += self.knownRxTypes[self.rxType][1]
P = self.getP(mesh)

if not adjoint:
Expand Down Expand Up @@ -185,9 +154,7 @@ def projectFields(self, u):
data = SimPEG.Survey.Data(self)
for src in self.srcList:
for rx in src.rxList:
print rx.nD
dat = rx.projectFields(src, self.mesh, u)
print dat.shape
data[src, rx] = rx.projectFields(src, self.mesh, u)
return data

Expand Down
13 changes: 13 additions & 0 deletions SimPEG/Mesh/TensorMesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -234,6 +234,9 @@ def getInterpolationMat(self, loc, locType='CC', zerosOutside=False):
'Fz' -> z-component of field defined on faces
'N' -> scalar field defined on nodes
'CC' -> scalar field defined on cell centers
'CCVx' -> x-component of vector field defined on cell centers
'CCVy' -> y-component of vector field defined on cell centers
'CCVz' -> z-component of vector field defined on cell centers
"""
if self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']:
raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType)
Expand All @@ -257,6 +260,16 @@ def getInterpolationMat(self, loc, locType='CC', zerosOutside=False):
Q = sp.hstack(components)
elif locType in ['CC', 'N']:
Q = Utils.interpmat(loc, *self.getTensor(locType))
elif locType in ['CCVx', 'CCVy', 'CCVz']:
Q = Utils.interpmat(loc, *self.getTensor('CC'))
Z = Utils.spzeros(loc.shape[0],self.nC)
if locType == 'CCVx':
Q = sp.hstack([Q,Z,Z])
elif locType == 'CCVy':
Q = sp.hstack([Z,Q,Z])
elif locType == 'CCVz':
Q = sp.hstack([Z,Z,Q])

else:
raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))

Expand Down
49 changes: 24 additions & 25 deletions tests/em/fdem/forward/test_FDEM_forwardEJHB.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@
testHJ = True
testEJ = True
testBH = True
verbose = False

TOLEBHJ = 1e-5
TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.)
Expand All @@ -21,56 +20,56 @@
class FDEM_CrossCheck(unittest.TestCase):
if testEJ:
def test_EJ_CrossCheck_jxr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxr', TOL=TOLEJHB))
def test_EJ_CrossCheck_jyr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyr', TOL=TOLEJHB))
def test_EJ_CrossCheck_jzr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzr', TOL=TOLEJHB))
def test_EJ_CrossCheck_jxi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxi', TOL=TOLEJHB))
def test_EJ_CrossCheck_jyi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyi', TOL=TOLEJHB))
def test_EJ_CrossCheck_jzi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzi', TOL=TOLEJHB))

def test_EJ_CrossCheck_exr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exr', TOL=TOLEJHB))
def test_EJ_CrossCheck_eyr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyr', TOL=TOLEJHB))
def test_EJ_CrossCheck_ezr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezr', TOL=TOLEJHB))
def test_EJ_CrossCheck_exi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exi', TOL=TOLEJHB))
def test_EJ_CrossCheck_eyi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyi', TOL=TOLEJHB))
def test_EJ_CrossCheck_ezi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezi', TOL=TOLEJHB))

def test_EJ_CrossCheck_bxr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxr', TOL=TOLEJHB))
def test_EJ_CrossCheck_byr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byr', TOL=TOLEJHB))
def test_EJ_CrossCheck_bzr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzr', TOL=TOLEJHB))
def test_EJ_CrossCheck_bxi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxi', TOL=TOLEJHB))
def test_EJ_CrossCheck_byi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byi', TOL=TOLEJHB))
def test_EJ_CrossCheck_bzi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzi', TOL=TOLEJHB))

def test_EJ_CrossCheck_hxr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxr', TOL=TOLEJHB))
def test_EJ_CrossCheck_hyr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyr', TOL=TOLEJHB))
def test_EJ_CrossCheck_hzr_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzr', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzr', TOL=TOLEJHB))
def test_EJ_CrossCheck_hxi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxi', TOL=TOLEJHB))
def test_EJ_CrossCheck_hyi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyi', TOL=TOLEJHB))
def test_EJ_CrossCheck_hzi_Jform(self):
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzi', verbose=verbose, TOL=TOLEJHB))
self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzi', TOL=TOLEJHB))

if __name__ == '__main__':
unittest.main()

0 comments on commit 73edef5

Please sign in to comment.