mirror of
https://github.com/wassname/simpeg.git
synced 2026-09-13 13:03:14 +08:00
Compare commits
17
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
8f1df12252 | ||
|
|
002554209b | ||
|
|
76f140f89b | ||
|
|
5ce26a97ae | ||
|
|
a087c6e5b1 | ||
|
|
44ec0623d6 | ||
|
|
57d0c84b41 | ||
|
|
a7befd6783 | ||
|
|
7acf46ddb8 | ||
|
|
b54d2494f5 | ||
|
|
45f5906554 | ||
|
|
8093288391 | ||
|
|
7457912d84 | ||
|
|
910d4214b3 | ||
|
|
20e80ed983 | ||
|
|
3e8580f28d | ||
|
|
3ad2c98d43 |
+8
-5
@@ -21,17 +21,20 @@ SimPEG
|
||||
:target: https://travis-ci.org/simpeg/simpeg
|
||||
:alt: Travis CI build status
|
||||
|
||||
.. image:: https://img.shields.io/coveralls/simpeg/simpeg.svg
|
||||
:target: https://coveralls.io/r/simpeg/simpeg?branch=master
|
||||
:alt: Coverage status
|
||||
|
||||
.. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
|
||||
:alt: gitter chat room at https://gitter.im/simpeg/simpeg
|
||||
:target: https://gitter.im/simpeg/simpeg
|
||||
|
||||
.. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg
|
||||
:target: https://codecov.io/gh/simpeg/simpeg
|
||||
:target: https://codecov.io/gh/simpeg/simpeg
|
||||
|
||||
.. image:: https://www.quantifiedcode.com/api/v1/project/933aa3decf444538aa432c8817169b6d/badge.svg
|
||||
:target: https://www.quantifiedcode.com/app/project/933aa3decf444538aa432c8817169b6d
|
||||
:alt: Code issues
|
||||
|
||||
.. image:: https://api.codacy.com/project/badge/Grade/4fc959a5294a418fa21fc7bc3b3aa078
|
||||
:target: https://www.codacy.com/app/lindseyheagy/simpeg?utm_source=github.com&utm_medium=referral&utm_content=simpeg/simpeg&utm_campaign=Badge_Grade
|
||||
:alt: codacy
|
||||
|
||||
Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
|
||||
|
||||
|
||||
@@ -476,7 +476,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
|
||||
fid.write('! ' + surveyType + ' FORMAT\n')
|
||||
|
||||
if iptype!=0:
|
||||
fid.write('IPTYPE=%i\n'%iptype)
|
||||
fid.write('IPTYPE={0:d}\n'.format(iptype))
|
||||
|
||||
else:
|
||||
fid.write('! ' + stype + ' FORMAT\n')
|
||||
@@ -512,7 +512,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
|
||||
|
||||
if surveyType == 'SURFACE':
|
||||
|
||||
fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
|
||||
fid.writelines("{0:f} ".format(ii) for ii in mkvc(tx[0,:]))
|
||||
M = M[:,0]
|
||||
N = N[:,0]
|
||||
|
||||
@@ -521,7 +521,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
|
||||
# Flip sign for z-elevation to depth
|
||||
tx[2::2,:] = -tx[2::2,:]
|
||||
|
||||
fid.writelines("%e " % ii for ii in mkvc(tx[::2,:]))
|
||||
fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[::2,:]))
|
||||
M = M[:,0::2]
|
||||
N = N[:,0::2]
|
||||
|
||||
@@ -529,22 +529,22 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
|
||||
M[:,1::2] = -M[:,1::2]
|
||||
N[:,1::2] = -N[:,1::2]
|
||||
|
||||
fid.write('%i\n'% nD)
|
||||
fid.write('{0:d}\n'.format(nD))
|
||||
np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
|
||||
|
||||
if dim=='3D':
|
||||
|
||||
if surveyType == 'SURFACE':
|
||||
|
||||
fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:]))
|
||||
fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[0:2,:]))
|
||||
M = M[:,0:2]
|
||||
N = N[:,0:2]
|
||||
|
||||
if surveyType == 'GENERAL':
|
||||
|
||||
fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
|
||||
fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[0:3,:]))
|
||||
|
||||
fid.write('%i\n'% nD)
|
||||
fid.write('{0:d}\n'.format(nD))
|
||||
np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
|
||||
fid.write('\n')
|
||||
|
||||
|
||||
+14
-14
@@ -15,7 +15,7 @@ class InversionDirective(object):
|
||||
@inversion.setter
|
||||
def inversion(self, i):
|
||||
if getattr(self,'_inversion',None) is not None:
|
||||
print 'Warning: InversionDirective %s has switched to a new inversion.' % self.__name__
|
||||
print 'Warning: InversionDirective {0!s} has switched to a new inversion.'.format(self.__name__)
|
||||
self._inversion = i
|
||||
|
||||
@property
|
||||
@@ -47,7 +47,7 @@ class DirectiveList(object):
|
||||
def __init__(self, *directives, **kwargs):
|
||||
self.dList = []
|
||||
for d in directives:
|
||||
assert isinstance(d, InversionDirective), 'All directives must be InversionDirectives not %s' % d.__name__
|
||||
assert isinstance(d, InversionDirective), 'All directives must be InversionDirectives not {0!s}'.format(d.__name__)
|
||||
self.dList.append(d)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
@@ -68,7 +68,7 @@ class DirectiveList(object):
|
||||
def inversion(self, i):
|
||||
if self.inversion is i: return
|
||||
if getattr(self,'_inversion',None) is not None:
|
||||
print 'Warning: %s has switched to a new inversion.' % self.__name__
|
||||
print 'Warning: {0!s} has switched to a new inversion.'.format(self.__name__)
|
||||
for d in self.dList:
|
||||
d.inversion = i
|
||||
self._inversion = i
|
||||
@@ -79,7 +79,7 @@ class DirectiveList(object):
|
||||
return
|
||||
|
||||
directives = ['initialize', 'endIter', 'finish']
|
||||
assert ruleType in directives, 'Directive type must be in ["%s"]' % '", "'.join(directives)
|
||||
assert ruleType in directives, 'Directive type must be in ["{0!s}"]'.format('", "'.join(directives))
|
||||
for r in self.dList:
|
||||
getattr(r, ruleType)()
|
||||
|
||||
@@ -141,7 +141,7 @@ class BetaSchedule(InversionDirective):
|
||||
|
||||
def endIter(self):
|
||||
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
|
||||
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
|
||||
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: {0:d}'.format(self.opt.iter)
|
||||
self.invProb.beta /= self.coolingFactor
|
||||
|
||||
|
||||
@@ -181,7 +181,7 @@ class SaveEveryIteration(InversionDirective):
|
||||
def fileName(self):
|
||||
if getattr(self, '_fileName', None) is None:
|
||||
from datetime import datetime
|
||||
self._fileName = '%s-%s'%(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M'))
|
||||
self._fileName = '{0!s}-{1!s}'.format(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M'))
|
||||
return self._fileName
|
||||
@fileName.setter
|
||||
def fileName(self, value):
|
||||
@@ -192,31 +192,31 @@ class SaveModelEveryIteration(SaveEveryIteration):
|
||||
"""SaveModelEveryIteration"""
|
||||
|
||||
def initialize(self):
|
||||
print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName
|
||||
print "SimPEG.SaveModelEveryIteration will save your models as: '###-{0!s}.npy'".format(self.fileName)
|
||||
|
||||
def endIter(self):
|
||||
np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
|
||||
np.save('{0:03d}-{1!s}'.format(self.opt.iter, self.fileName), self.opt.xc)
|
||||
|
||||
|
||||
class SaveOutputEveryIteration(SaveEveryIteration):
|
||||
"""SaveModelEveryIteration"""
|
||||
|
||||
def initialize(self):
|
||||
print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName
|
||||
print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-{0!s}.txt'".format(self.fileName)
|
||||
f = open(self.fileName+'.txt', 'w')
|
||||
f.write(" # beta phi_d phi_m f\n")
|
||||
f.close()
|
||||
|
||||
def endIter(self):
|
||||
f = open(self.fileName+'.txt', 'a')
|
||||
f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
|
||||
f.write(' {0:3d} {1:1.4e} {2:1.4e} {3:1.4e} {4:1.4e}\n'.format(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
|
||||
f.close()
|
||||
|
||||
class SaveOutputDictEveryIteration(SaveEveryIteration):
|
||||
"""SaveOutputDictEveryIteration"""
|
||||
|
||||
def initialize(self):
|
||||
print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
|
||||
print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-{0!s}.npz'".format(self.fileName)
|
||||
|
||||
def endIter(self):
|
||||
# Save the data.
|
||||
@@ -328,7 +328,7 @@ class Update_IRLS(InversionDirective):
|
||||
|
||||
# Beta Schedule
|
||||
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
|
||||
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
|
||||
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: {0:d}'.format(self.opt.iter)
|
||||
self.invProb.beta /= self.coolingFactor
|
||||
|
||||
|
||||
@@ -340,11 +340,11 @@ class Update_IRLS(InversionDirective):
|
||||
phim_new = self.reg.eval(self.invProb.curModel)
|
||||
self.f_change = np.abs(self.f_old - phim_new) / self.f_old
|
||||
|
||||
print "Regularization decrease: %6.3e" % (self.f_change)
|
||||
print "Regularization decrease: {0:6.3e}".format((self.f_change))
|
||||
|
||||
# Check for maximum number of IRLS cycles
|
||||
if self.IRLSiter == self.maxIRLSiter:
|
||||
print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
|
||||
print "Reach maximum number of IRLS cycles: {0:d}".format(self.maxIRLSiter)
|
||||
self.opt.stopNextIteration = True
|
||||
return
|
||||
|
||||
|
||||
@@ -42,7 +42,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: total electric field
|
||||
"""
|
||||
if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
|
||||
raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting e from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
|
||||
|
||||
@@ -56,7 +56,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: total magnetic flux density
|
||||
"""
|
||||
if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
|
||||
raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting b from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
|
||||
|
||||
@@ -70,7 +70,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: total magnetic field
|
||||
"""
|
||||
if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
|
||||
raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting h from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
|
||||
|
||||
@@ -84,7 +84,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: total current density
|
||||
"""
|
||||
if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
|
||||
raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting j from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
|
||||
|
||||
@@ -100,7 +100,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: derivative times a vector (or tuple for adjoint)
|
||||
"""
|
||||
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
if adjoint:
|
||||
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
|
||||
@@ -118,7 +118,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: derivative times a vector (or tuple for adjoint)
|
||||
"""
|
||||
if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting bDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
if adjoint:
|
||||
return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
|
||||
@@ -136,7 +136,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: derivative times a vector (or tuple for adjoint)
|
||||
"""
|
||||
if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting hDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
if adjoint:
|
||||
return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
|
||||
@@ -154,7 +154,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: derivative times a vector (or tuple for adjoint)
|
||||
"""
|
||||
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
if adjoint:
|
||||
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
|
||||
|
||||
@@ -11,8 +11,8 @@ class BaseRx(SimPEG.Survey.BaseRx):
|
||||
"""
|
||||
|
||||
def __init__(self, locs, orientation=None, component=None):
|
||||
assert(orientation in ['x','y','z']), "Orientation %s not known. Orientation must be in 'x', 'y', 'z'. Arbitrary orientations have not yet been implemented."%orientation
|
||||
assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not %s"%component
|
||||
assert(orientation in ['x','y','z']), "Orientation {0!s} not known. Orientation must be in 'x', 'y', 'z'. Arbitrary orientations have not yet been implemented.".format(orientation)
|
||||
assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not {0!s}".format(component)
|
||||
|
||||
self.projComp = orientation
|
||||
self.component = component
|
||||
|
||||
@@ -9,7 +9,7 @@ class Fields(SimPEG.Problem.Fields):
|
||||
|
||||
def _phiDeriv(self, src, du_dm_v, v, adjoint=False):
|
||||
if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting phiDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
if adjoint:
|
||||
return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
|
||||
@@ -18,7 +18,7 @@ class Fields(SimPEG.Problem.Fields):
|
||||
|
||||
def _eDeriv(self, src, du_dm_v, v, adjoint=False):
|
||||
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
if adjoint:
|
||||
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
|
||||
@@ -26,7 +26,7 @@ class Fields(SimPEG.Problem.Fields):
|
||||
|
||||
def _jDeriv(self, src, du_dm_v, v, adjoint=False):
|
||||
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
if adjoint:
|
||||
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
|
||||
|
||||
@@ -32,7 +32,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
|
||||
|
||||
def _phiDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
|
||||
if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting phiDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
if adjoint:
|
||||
return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
|
||||
@@ -41,7 +41,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
|
||||
|
||||
def _eDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
|
||||
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
if adjoint:
|
||||
return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
|
||||
@@ -49,7 +49,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
|
||||
|
||||
def _jDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
|
||||
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
|
||||
|
||||
if adjoint:
|
||||
return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)
|
||||
|
||||
@@ -46,7 +46,7 @@ class BaseDCProblem(BaseEMProblem):
|
||||
du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
|
||||
|
||||
for rx in src.rxList:
|
||||
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
|
||||
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
|
||||
Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
|
||||
return Utils.mkvc(Jv)
|
||||
@@ -69,7 +69,7 @@ class BaseDCProblem(BaseEMProblem):
|
||||
u_src = f[src, self._solutionType]
|
||||
for rx in src.rxList:
|
||||
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
|
||||
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
|
||||
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
|
||||
|
||||
ATinvdf_duT = self.Ainv * df_duT
|
||||
|
||||
@@ -60,7 +60,7 @@ class BaseDCProblem_2D(BaseEMProblem):
|
||||
dRHS_dm_v = self.getRHSDeriv(ky, src, v)
|
||||
du_dm_v = self.Ainv[iky] * ( - dA_dm_v + dRHS_dm_v )
|
||||
for rx in src.rxList:
|
||||
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
|
||||
df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False)
|
||||
# Trapezoidal intergration
|
||||
Jv1_temp = 1./np.pi*rx.evalDeriv(ky, src, self.mesh, f, df_dm_v)
|
||||
@@ -101,7 +101,7 @@ class BaseDCProblem_2D(BaseEMProblem):
|
||||
ky = self.kys[iky]
|
||||
AT = self.getA(ky)
|
||||
PTv = rx.evalDeriv(ky, src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
|
||||
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
|
||||
df_duT, df_dmT = df_duTFun(iky, src, None, PTv, adjoint=True)
|
||||
|
||||
ATinvdf_duT = self.Ainv[iky] * df_duT
|
||||
|
||||
@@ -56,7 +56,7 @@ class BaseIPProblem(BaseEMProblem):
|
||||
du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
|
||||
|
||||
for rx in src.rxList:
|
||||
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
|
||||
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
|
||||
Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
|
||||
# Conductivity (d u / d log sigma)
|
||||
@@ -83,7 +83,7 @@ class BaseIPProblem(BaseEMProblem):
|
||||
u_src = f[src, self._solutionType]
|
||||
for rx in src.rxList:
|
||||
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
|
||||
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
|
||||
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
|
||||
ATinvdf_duT = self.Ainv * df_duT
|
||||
dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
|
||||
|
||||
@@ -83,7 +83,7 @@ class BaseSIPProblem(BaseEMProblem):
|
||||
for rx in src.rxList:
|
||||
timeindex = rx.getTimeP(self.survey.times)
|
||||
if timeindex[tind]:
|
||||
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
|
||||
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
|
||||
Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
|
||||
|
||||
@@ -122,7 +122,7 @@ class BaseSIPProblem(BaseEMProblem):
|
||||
for rx in src.rxList:
|
||||
timeindex = rx.getTimeP(self.survey.times)
|
||||
if timeindex[tind]:
|
||||
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
|
||||
df_dm_v0 = df_dmFun(src, du_dm_v0, v0, adjoint=False)
|
||||
df_dm_v1 = df_dmFun(src, du_dm_v1, v1, adjoint=False)
|
||||
Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v0)
|
||||
@@ -153,7 +153,7 @@ class BaseSIPProblem(BaseEMProblem):
|
||||
timeindex = rx.getTimeP(self.survey.times)
|
||||
if timeindex[tind]:
|
||||
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx, t], adjoint=True) # wrt f, need possibility wrt m
|
||||
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
|
||||
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
|
||||
ATinvdf_duT = self.Ainv * df_duT
|
||||
dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
|
||||
|
||||
+10
-10
@@ -47,7 +47,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
self.waveformType = "GENERAL"
|
||||
|
||||
def fields(self, m):
|
||||
if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print '{0!s}\nCalculating fields(m)\n{1!s}'.format('*'*50, '*'*50)
|
||||
self.curModel = m
|
||||
# Create a fields storage object
|
||||
F = self._FieldsForward_pair(self.mesh, self.survey)
|
||||
@@ -55,7 +55,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
# Set the initial conditions
|
||||
F[src,:,0] = src.getInitialFields(self.mesh)
|
||||
F = self.forward(m, self.getRHS, F=F)
|
||||
if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print '{0!s}\nDone calculating fields(m)\n{1!s}'.format('*'*50, '*'*50)
|
||||
return F
|
||||
|
||||
def forward(self, m, RHS, F=None):
|
||||
@@ -70,11 +70,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
if Ainv is not None:
|
||||
Ainv.clean()
|
||||
A = self.getA(tInd)
|
||||
if self.verbose: print 'Factoring... (dt = %e)'%dt
|
||||
if self.verbose: print 'Factoring... (dt = {0:e})'.format(dt)
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
if self.verbose: print 'Done'
|
||||
rhs = RHS(tInd, F)
|
||||
if self.verbose: print ' Solving... (tInd = %d)'%tInd
|
||||
if self.verbose: print ' Solving... (tInd = {0:d})'.format(tInd)
|
||||
sol = Ainv * rhs
|
||||
if self.verbose: print ' Done...'
|
||||
if sol.ndim == 1:
|
||||
@@ -95,11 +95,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
if Ainv is not None:
|
||||
Ainv.clean()
|
||||
A = self.getA(tInd)
|
||||
if self.verbose: print 'Factoring (Adjoint)... (dt = %e)'%dt
|
||||
if self.verbose: print 'Factoring (Adjoint)... (dt = {0:e})'.format(dt)
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
if self.verbose: print 'Done'
|
||||
rhs = RHS(tInd, F)
|
||||
if self.verbose: print ' Solving (Adjoint)... (tInd = %d)'%tInd
|
||||
if self.verbose: print ' Solving (Adjoint)... (tInd = {0:d})'.format(tInd)
|
||||
sol = Ainv * rhs
|
||||
if self.verbose: print ' Done...'
|
||||
if sol.ndim == 1:
|
||||
@@ -123,14 +123,14 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
* Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
|
||||
|
||||
"""
|
||||
if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print '{0!s}\nCalculating J(v)\n{1!s}'.format('*'*50, '*'*50)
|
||||
self.curModel = m
|
||||
if f is None:
|
||||
f = self.fields(m)
|
||||
p = self.Gvec(m, v, f)
|
||||
y = self.solveAh(m, p)
|
||||
Jv = self.survey.evalDeriv(f, v=y)
|
||||
if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print '{0!s}\nDone calculating J(v)\n{1!s}'.format('*'*50, '*'*50)
|
||||
return - mkvc(Jv)
|
||||
|
||||
def Jtvec(self, m, v, f=None):
|
||||
@@ -148,7 +148,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
* Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\)
|
||||
|
||||
"""
|
||||
if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print '{0!s}\nCalculating J^T(v)\n{1!s}'.format('*'*50, '*'*50)
|
||||
self.curModel = m
|
||||
if f is None:
|
||||
f = self.fields(m)
|
||||
@@ -159,6 +159,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
p = self.survey.evalDeriv(f, v=v, adjoint=True)
|
||||
y = self.solveAht(m, p)
|
||||
w = self.Gtvec(m, y, f)
|
||||
if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print '{0!s}\nDone calculating J^T(v)\n{1!s}'.format('*'*50, '*'*50)
|
||||
return - mkvc(w)
|
||||
|
||||
|
||||
@@ -58,7 +58,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
|
||||
Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*S_m, S_e))
|
||||
|
||||
if verbose:
|
||||
print ' Fetching %s problem' % (fdemType)
|
||||
print ' Fetching {0!s} problem'.format((fdemType))
|
||||
|
||||
if fdemType == 'e':
|
||||
survey = EM.FDEM.Survey(Src)
|
||||
@@ -94,7 +94,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
|
||||
|
||||
prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
|
||||
mesh = prb1.mesh
|
||||
print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
|
||||
print 'Cross Checking Forward: {0!s}, {1!s} formulations - {2!s}'.format(fdemType1, fdemType2, comp)
|
||||
|
||||
logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
|
||||
mu = np.ones(mesh.nC)*MU
|
||||
|
||||
@@ -110,7 +110,7 @@ def run(plotIt=True):
|
||||
# Mesh
|
||||
mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
|
||||
|
||||
print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
|
||||
print 'Mesh Extent xmax: {0:f},: zmin: {1:f}, zmax: {2:f}'.format(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
|
||||
print 'Number of cells', mesh.nC
|
||||
|
||||
if plotIt is True:
|
||||
|
||||
@@ -98,7 +98,7 @@ def run(plotIt=True, n=60):
|
||||
ii = int(ii)
|
||||
out = M.plotImage(PHIS[ii][1],ax=ax)
|
||||
ax.axis('off')
|
||||
ax.set_title('Elapsed Time: %4.1f'%PHIS[ii][0])
|
||||
ax.set_title('Elapsed Time: {0:4.1f}'.format(PHIS[ii][0]))
|
||||
plt.show()
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
@@ -29,15 +29,15 @@ def run(plotIt=True, n=60):
|
||||
axes[0].set_ylim([-1,17])
|
||||
|
||||
for ii, loc in zip(range(M.nC),M.gridCC):
|
||||
axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='r')
|
||||
axes[0].text(loc[0]+0.2,loc[1],'{0:d}'.format(ii), color='r')
|
||||
|
||||
axes[0].plot(M.gridFx[:,0],M.gridFx[:,1], 'g>')
|
||||
for ii, loc in zip(range(M.nFx),M.gridFx):
|
||||
axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='g')
|
||||
axes[0].text(loc[0]+0.2,loc[1],'{0:d}'.format(ii), color='g')
|
||||
|
||||
axes[0].plot(M.gridFy[:,0],M.gridFy[:,1], 'm^')
|
||||
for ii, loc in zip(range(M.nFy),M.gridFy):
|
||||
axes[0].text(loc[0]+0.2,loc[1]+0.2,'%d'%(ii+M.nFx), color='m')
|
||||
axes[0].text(loc[0]+0.2,loc[1]+0.2,'{0:d}'.format((ii+M.nFx)), color='m')
|
||||
|
||||
axes[1].spy(M.faceDiv)
|
||||
axes[1].set_title('Face Divergence')
|
||||
|
||||
@@ -59,7 +59,7 @@ if __name__ == '__main__':
|
||||
if line == "##### AUTOIMPORTS #####\n":
|
||||
inimports = not inimports
|
||||
if inimports:
|
||||
out += '\n'.join(["import %s"%_ for _ in exfiles])
|
||||
out += '\n'.join(["import {0!s}".format(_) for _ in exfiles])
|
||||
out += '\n\n__examples__ = ["' + '", "'.join(exfiles)+ '"]\n'
|
||||
out += '\n##### AUTOIMPORTS #####\n'
|
||||
f.close()
|
||||
@@ -76,11 +76,11 @@ if __name__ == '__main__':
|
||||
|
||||
docstr = runFunction.__doc__
|
||||
if docstr is None:
|
||||
doc = '%s\n%s'%(name.replace('_',' '),'='*len(name))
|
||||
doc = '{0!s}\n{1!s}'.format(name.replace('_',' '), '='*len(name))
|
||||
else:
|
||||
doc = '\n'.join([_[8:].rstrip() for _ in docstr.split('\n')])
|
||||
|
||||
out = """.. _examples_%s:
|
||||
out = """.. _examples_{0!s}:
|
||||
|
||||
.. --------------------------------- ..
|
||||
.. ..
|
||||
@@ -90,21 +90,21 @@ if __name__ == '__main__':
|
||||
.. ..
|
||||
.. --------------------------------- ..
|
||||
|
||||
%s
|
||||
{1!s}
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.%s.run()
|
||||
Examples.{2!s}.run()
|
||||
|
||||
.. literalinclude:: ../../../SimPEG/Examples/%s.py
|
||||
.. literalinclude:: ../../../SimPEG/Examples/{3!s}.py
|
||||
:language: python
|
||||
:linenos:
|
||||
"""%(name,doc,name,name)
|
||||
""".format(name, doc, name, name)
|
||||
|
||||
rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst']))
|
||||
|
||||
print 'Creating: %s.rst'%name
|
||||
print 'Creating: {0!s}.rst'.format(name)
|
||||
f = open(rst, 'w')
|
||||
f.write(out)
|
||||
f.close()
|
||||
|
||||
@@ -116,7 +116,7 @@ class RichardsMap(object):
|
||||
ax.semilogx(self.k(h, m), h)
|
||||
|
||||
def _assertMatchesPair(self, pair):
|
||||
assert isinstance(self, pair), "Mapping object must be an instance of a %s class."%(pair.__name__)
|
||||
assert isinstance(self, pair), "Mapping object must be an instance of a {0!s} class.".format((pair.__name__))
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -140,7 +140,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
|
||||
for ii, dt in enumerate(self.timeSteps):
|
||||
bc = self.getBoundaryConditions(ii, u[ii])
|
||||
u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, dt, bc, return_g=return_g), u[ii])
|
||||
if self.debug: print "Solving Fields (%4d/%d - %3.1f%% Done) %d Iterations, %4.2f seconds"%(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic)
|
||||
if self.debug: print "Solving Fields ({0:4d}/{1:d} - {2:3.1f}% Done) {3:d} Iterations, {4:4.2f} seconds".format(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic)
|
||||
return u
|
||||
|
||||
@Utils.timeIt
|
||||
|
||||
+4
-4
@@ -37,7 +37,7 @@ class Fields(object):
|
||||
for f in self.knownFields:
|
||||
loc =self.knownFields[f]
|
||||
sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2)
|
||||
return "%e MB"%sz
|
||||
return "{0:e} MB".format(sz)
|
||||
|
||||
def _storageShape(self, loc):
|
||||
nSrc = self.survey.nSrc
|
||||
@@ -84,12 +84,12 @@ class Fields(object):
|
||||
return
|
||||
if accessType=='set' and name not in self.knownFields:
|
||||
if name in self.aliasFields:
|
||||
raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name)
|
||||
raise KeyError("Invalid field name ({0!s}) for setter, you can't set an aliased property".format(name))
|
||||
else:
|
||||
raise KeyError('Invalid field name (%s) for setter'%name)
|
||||
raise KeyError('Invalid field name ({0!s}) for setter'.format(name))
|
||||
|
||||
elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields):
|
||||
raise KeyError('Invalid field name (%s) for getter'%name)
|
||||
raise KeyError('Invalid field name ({0!s}) for getter'.format(name))
|
||||
return name
|
||||
|
||||
def _indexAndNameFromKey(self, key, accessType):
|
||||
|
||||
+7
-7
@@ -101,7 +101,7 @@ class IdentityMap(object):
|
||||
:return: passed the test?
|
||||
|
||||
"""
|
||||
print 'Testing %s' % str(self)
|
||||
print 'Testing {0!s}'.format(str(self))
|
||||
if m is None:
|
||||
m = abs(np.random.rand(self.nP))
|
||||
if 'plotIt' not in kwargs:
|
||||
@@ -111,21 +111,21 @@ class IdentityMap(object):
|
||||
def _assertMatchesPair(self, pair):
|
||||
assert (isinstance(self, pair) or
|
||||
isinstance(self, ComboMap) and isinstance(self.maps[0], pair)
|
||||
), "Mapping object must be an instance of a %s class."%(pair.__name__)
|
||||
), "Mapping object must be an instance of a {0!s} class.".format((pair.__name__))
|
||||
|
||||
def __mul__(self, val):
|
||||
if isinstance(val, IdentityMap):
|
||||
if not (self.shape[1] == '*' or val.shape[0] == '*') and not self.shape[1] == val.shape[0]:
|
||||
raise ValueError('Dimension mismatch in %s and %s.' % (str(self), str(val)))
|
||||
raise ValueError('Dimension mismatch in {0!s} and {1!s}.'.format(str(self), str(val)))
|
||||
return ComboMap([self, val])
|
||||
elif isinstance(val, np.ndarray):
|
||||
if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]:
|
||||
raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape)))
|
||||
raise ValueError('Dimension mismatch in {0!s} and np.ndarray{1!s}.'.format(str(self), str(val.shape)))
|
||||
return self._transform(val)
|
||||
raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!')
|
||||
|
||||
def __str__(self):
|
||||
return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
|
||||
return "{0!s}({1!s},{2!s})".format(self.__class__.__name__, self.shape[0], self.shape[1])
|
||||
|
||||
|
||||
class ComboMap(IdentityMap):
|
||||
@@ -140,7 +140,7 @@ class ComboMap(IdentityMap):
|
||||
if ii > 0 and not (self.shape[1] == '*' or m.shape[0] == '*') and not self.shape[1] == m.shape[0]:
|
||||
prev = self.maps[-1]
|
||||
errArgs = (prev.__class__.__name__, prev.shape[0], prev.shape[1], m.__class__.__name__, m.shape[0], m.shape[1])
|
||||
raise ValueError('Dimension mismatch in map[%s] (%s, %s) and map[%s] (%s, %s).' % errArgs)
|
||||
raise ValueError('Dimension mismatch in map[{0!s}] ({1!s}, {2!s}) and map[{3!s}] ({4!s}, {5!s}).'.format(*errArgs))
|
||||
|
||||
if isinstance(m, ComboMap):
|
||||
self.maps += m.maps
|
||||
@@ -173,7 +173,7 @@ class ComboMap(IdentityMap):
|
||||
return deriv
|
||||
|
||||
def __str__(self):
|
||||
return 'ComboMap[%s](%s,%s)' % (' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1])
|
||||
return 'ComboMap[{0!s}]({1!s},{2!s})'.format(' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1])
|
||||
|
||||
|
||||
class ExpMap(IdentityMap):
|
||||
|
||||
@@ -522,7 +522,7 @@ class BaseRectangularMesh(BaseMesh):
|
||||
assert xType in outType, 'You cannot change type of components.'
|
||||
if type(x) == list:
|
||||
for i, xi in enumerate(x):
|
||||
assert isinstance(x, np.ndarray), "x[%i] must be a numpy array" % i
|
||||
assert isinstance(x, np.ndarray), "x[{0:d}] must be a numpy array".format(i)
|
||||
assert xi.size == x[0].size, "Number of elements in list must not change."
|
||||
|
||||
x_array = np.ones((x.size, len(x)))
|
||||
|
||||
+166
-134
@@ -4,14 +4,28 @@ from DiffOperators import DiffOperators
|
||||
from InnerProducts import InnerProducts
|
||||
from View import CurvView
|
||||
|
||||
|
||||
# Some helper functions.
|
||||
length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
|
||||
length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
|
||||
normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
|
||||
normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
|
||||
def length2D(x):
|
||||
return (x[:, 0]**2 + x[:, 1]**2)**0.5
|
||||
|
||||
|
||||
class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView):
|
||||
def length3D(x):
|
||||
return (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
|
||||
|
||||
|
||||
def normalize2D(x):
|
||||
return x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
|
||||
|
||||
|
||||
def normalize3D(x):
|
||||
return x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
|
||||
|
||||
|
||||
# Curvi Mesh
|
||||
|
||||
class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts,
|
||||
CurvView):
|
||||
"""
|
||||
CurvilinearMesh is a mesh class that deals with curvilinear meshes.
|
||||
|
||||
@@ -31,12 +45,16 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
|
||||
_meshType = 'Curv'
|
||||
|
||||
def __init__(self, nodes):
|
||||
assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray"
|
||||
assert type(nodes) == list, ("'nodes' variable must be a list of "
|
||||
"np.ndarray")
|
||||
assert len(nodes) > 1, "len(node) must be greater than 1"
|
||||
|
||||
for i, nodes_i in enumerate(nodes):
|
||||
assert isinstance(nodes_i, np.ndarray), ("nodes[%i] is not a numpy array." % i)
|
||||
assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % i)
|
||||
assert isinstance(nodes_i, np.ndarray), ("nodes[{0:d}] is not a"
|
||||
"numpy array.".format(i))
|
||||
assert nodes_i.shape == nodes[0].shape, ("nodes[{0:d}] is not the "
|
||||
"same shape as nodes[0]"
|
||||
.format(i))
|
||||
|
||||
assert len(nodes[0].shape) == len(nodes), "Dimension mismatch"
|
||||
assert len(nodes[0].shape) > 1, "Not worth using Curv for a 1D mesh."
|
||||
@@ -48,80 +66,79 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
|
||||
for i, node_i in enumerate(nodes):
|
||||
self._gridN[:, i] = Utils.mkvc(node_i.astype(float))
|
||||
|
||||
def gridCC():
|
||||
doc = "Cell-centered grid."
|
||||
|
||||
def fget(self):
|
||||
if self._gridCC is None:
|
||||
self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:,i] for i in range(self.dim)]).reshape((-1,self.dim), order='F')
|
||||
@property
|
||||
def gridCC(self):
|
||||
"""
|
||||
Cell-centered grid
|
||||
"""
|
||||
if getattr(self, '_gridCC', None) is None:
|
||||
self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:, i]
|
||||
for i in range(self.dim)]).reshape(
|
||||
(-1, self.dim), order='F')
|
||||
return self._gridCC
|
||||
return locals()
|
||||
_gridCC = None # Store grid by default
|
||||
gridCC = property(**gridCC())
|
||||
|
||||
def gridN():
|
||||
doc = "Nodal grid."
|
||||
|
||||
def fget(self):
|
||||
if self._gridN is None:
|
||||
@property
|
||||
def gridN(self):
|
||||
"""
|
||||
Nodal grid.
|
||||
"""
|
||||
if getattr(self, '_gridN', None) is None:
|
||||
raise Exception("Someone deleted this. I blame you.")
|
||||
return self._gridN
|
||||
return locals()
|
||||
_gridN = None # Store grid by default
|
||||
gridN = property(**gridN())
|
||||
|
||||
def gridFx():
|
||||
doc = "Face staggered grid in the x direction."
|
||||
@property
|
||||
def gridFx(self):
|
||||
"""
|
||||
Face staggered grid in the x direction.
|
||||
"""
|
||||
|
||||
def fget(self):
|
||||
if self._gridFx is None:
|
||||
if getattr(self, '_gridFx', None) is None:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
|
||||
self._gridFx = np.c_[XY[0], XY[1]]
|
||||
elif self.dim == 3:
|
||||
XYZ = [Utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N]
|
||||
XYZ = [Utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] +
|
||||
n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N]
|
||||
self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
return self._gridFx
|
||||
return locals()
|
||||
_gridFx = None # Store grid by default
|
||||
gridFx = property(**gridFx())
|
||||
|
||||
def gridFy():
|
||||
doc = "Face staggered grid in the y direction."
|
||||
@property
|
||||
def gridFy(self):
|
||||
"""
|
||||
Face staggered grid in the y direction.
|
||||
"""
|
||||
|
||||
def fget(self):
|
||||
if self._gridFy is None:
|
||||
if getattr(self, '_gridFy', None) is None:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
|
||||
self._gridFy = np.c_[XY[0], XY[1]]
|
||||
elif self.dim == 3:
|
||||
XYZ = [Utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N]
|
||||
XYZ = [Utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] +
|
||||
n[1:, :, :-1] + n[1:, :, 1:])) for n in N]
|
||||
self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
return self._gridFy
|
||||
return locals()
|
||||
_gridFy = None # Store grid by default
|
||||
gridFy = property(**gridFy())
|
||||
|
||||
def gridFz():
|
||||
doc = "Face staggered grid in the z direction."
|
||||
@property
|
||||
def gridFz(self):
|
||||
"""
|
||||
Face staggered grid in the y direction.
|
||||
"""
|
||||
|
||||
def fget(self):
|
||||
if self._gridFz is None and self.dim == 3:
|
||||
if getattr(self, '_gridFz', None) is None:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
|
||||
XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] +
|
||||
n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
|
||||
self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
return self._gridFz
|
||||
return locals()
|
||||
_gridFz = None # Store grid by default
|
||||
gridFz = property(**gridFz())
|
||||
|
||||
def gridEx():
|
||||
doc = "Edge staggered grid in the x direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridEx is None:
|
||||
@property
|
||||
def gridEx(self):
|
||||
"""
|
||||
Edge staggered grid in the x direction.
|
||||
"""
|
||||
if getattr(self, '_gridEx', None) is None:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
|
||||
@@ -130,15 +147,13 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
|
||||
XYZ = [Utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N]
|
||||
self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
return self._gridEx
|
||||
return locals()
|
||||
_gridEx = None # Store grid by default
|
||||
gridEx = property(**gridEx())
|
||||
|
||||
def gridEy():
|
||||
doc = "Edge staggered grid in the y direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridEy is None:
|
||||
@property
|
||||
def gridEy(self):
|
||||
"""
|
||||
Edge staggered grid in the y direction.
|
||||
"""
|
||||
if getattr(self, '_gridEy', None) is None:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
if self.dim == 2:
|
||||
XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
|
||||
@@ -147,22 +162,17 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
|
||||
XYZ = [Utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N]
|
||||
self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
return self._gridEy
|
||||
return locals()
|
||||
_gridEy = None # Store grid by default
|
||||
gridEy = property(**gridEy())
|
||||
|
||||
def gridEz():
|
||||
doc = "Edge staggered grid in the z direction."
|
||||
|
||||
def fget(self):
|
||||
if self._gridEz is None and self.dim == 3:
|
||||
@property
|
||||
def gridEz(self):
|
||||
"""
|
||||
Edge staggered grid in the z direction.
|
||||
"""
|
||||
if getattr(self, '_gridEz', None) is None and self.dim == 3:
|
||||
N = self.r(self.gridN, 'N', 'N', 'M')
|
||||
XYZ = [Utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N]
|
||||
self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
|
||||
return self._gridEz
|
||||
return locals()
|
||||
_gridEz = None # Store grid by default
|
||||
gridEz = property(**gridEz())
|
||||
|
||||
# --------------- Geometries ---------------------
|
||||
#
|
||||
@@ -194,19 +204,25 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
|
||||
# | / | /
|
||||
# D -------------- C
|
||||
# node(i+1,j,k) node(i+1,j+1,k)
|
||||
def vol():
|
||||
doc = "Construct cell volumes of the 3D model as 1d array."
|
||||
|
||||
def fget(self):
|
||||
if(self._vol is None):
|
||||
@property
|
||||
def vol(self):
|
||||
"""
|
||||
Construct cell volumes of the 3D model as 1d array
|
||||
"""
|
||||
|
||||
if getattr(self, '_vol', None) is None:
|
||||
if self.dim == 2:
|
||||
A, B, C, D = Utils.indexCube('ABCD', self.vnC+1)
|
||||
normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
|
||||
normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros(
|
||||
(self.nN, 1))], A, B, C, D)
|
||||
self._vol = area
|
||||
elif self.dim == 3:
|
||||
# Each polyhedron can be decomposed into 5 tetrahedrons
|
||||
# However, this presents a choice so we may as well divide in two ways and average.
|
||||
A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.vnC+1)
|
||||
# However, this presents a choice so we may as well divide in
|
||||
# two ways and average.
|
||||
A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.vnC +
|
||||
1)
|
||||
|
||||
vol1 = (Utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top
|
||||
Utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top
|
||||
@@ -222,50 +238,60 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
|
||||
|
||||
self._vol = (vol1 + vol2)/2
|
||||
return self._vol
|
||||
return locals()
|
||||
_vol = None
|
||||
vol = property(**vol())
|
||||
|
||||
def area():
|
||||
doc = "Face areas."
|
||||
|
||||
def fget(self):
|
||||
if(self._area is None or self._normals is None):
|
||||
@property
|
||||
def area(self):
|
||||
if (getattr(self, '_area', None) is None or
|
||||
getattr(self, '_normals', None) is None):
|
||||
# Compute areas of cell faces
|
||||
if(self.dim == 2):
|
||||
xy = self.gridN
|
||||
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy]))
|
||||
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx,
|
||||
self.nCy]))
|
||||
edge1 = xy[B, :] - xy[A, :]
|
||||
normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
|
||||
area1 = length2D(edge1)
|
||||
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy]))
|
||||
# Note that we are doing A-D to make sure the normal points the right way.
|
||||
# Think about it. Look at the picture. Normal points towards C iff you do this.
|
||||
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
|
||||
self.nNy]))
|
||||
# Note that we are doing A-D to make sure the normal points the
|
||||
# right way.
|
||||
# Think about it. Look at the picture. Normal points towards C
|
||||
# iff you do this.
|
||||
edge2 = xy[A, :] - xy[D, :]
|
||||
normal2 = np.c_[edge2[:, 1], -edge2[:, 0]]
|
||||
area2 = length2D(edge2)
|
||||
self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)]
|
||||
self._normals = [normalize2D(normal1), normalize2D(normal2)]
|
||||
|
||||
elif(self.dim == 3):
|
||||
|
||||
A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array([self.nNx, self.nCy, self.nCz]))
|
||||
normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False)
|
||||
A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array(
|
||||
[self.nNx, self.nCy, self.nCz]))
|
||||
normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B,
|
||||
average=False,
|
||||
normalizeNormals=False)
|
||||
|
||||
A, D, H, E = Utils.indexCube('ADHE', self.vnC+1, np.array([self.nCx, self.nNy, self.nCz]))
|
||||
normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False)
|
||||
A, D, H, E = Utils.indexCube('ADHE', self.vnC+1, np.array(
|
||||
[self.nCx, self.nNy, self.nCz]))
|
||||
normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E,
|
||||
average=False,
|
||||
normalizeNormals=False)
|
||||
|
||||
A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array([self.nCx, self.nCy, self.nNz]))
|
||||
normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False)
|
||||
A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array(
|
||||
[self.nCx, self.nCy, self.nNz]))
|
||||
normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D,
|
||||
average=False,
|
||||
normalizeNormals=False)
|
||||
|
||||
self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)]
|
||||
self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2),
|
||||
Utils.mkvc(area3)]
|
||||
self._normals = [normal1, normal2, normal3]
|
||||
return self._area
|
||||
return locals()
|
||||
_area = None
|
||||
area = property(**area())
|
||||
|
||||
def normals():
|
||||
doc = """Face normals: calling this will average
|
||||
@property
|
||||
def normals(self):
|
||||
"""
|
||||
Face normals: calling this will average
|
||||
the computed normals so that there is one
|
||||
per face. This is especially relevant in
|
||||
3D, as there are up to 4 different normals
|
||||
@@ -276,8 +302,7 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
|
||||
NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M')
|
||||
"""
|
||||
|
||||
def fget(self):
|
||||
if(self._normals is None):
|
||||
if getattr(self, '_normals', None) is None:
|
||||
self.area # calling .area will create the face normals
|
||||
if self.dim == 2:
|
||||
return normalize2D(np.r_[self._normals[0], self._normals[1]])
|
||||
@@ -286,48 +311,55 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
|
||||
normal2 = (self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3])/4
|
||||
normal3 = (self._normals[2][0] + self._normals[2][1] + self._normals[2][2] + self._normals[2][3])/4
|
||||
return normalize3D(np.r_[normal1, normal2, normal3])
|
||||
return locals()
|
||||
_normals = None
|
||||
normals = property(**normals())
|
||||
|
||||
def edge():
|
||||
doc = "Edge legnths."
|
||||
|
||||
def fget(self):
|
||||
if(self._edge is None or self._tangents is None):
|
||||
@property
|
||||
def edge(self):
|
||||
"""
|
||||
Edge lengths
|
||||
"""
|
||||
if getattr(self, '_edge', None) is None:
|
||||
if(self.dim == 2):
|
||||
xy = self.gridN
|
||||
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy]))
|
||||
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
|
||||
self.nNy]))
|
||||
edge1 = xy[D, :] - xy[A, :]
|
||||
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy]))
|
||||
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx,
|
||||
self.nCy]))
|
||||
edge2 = xy[B, :] - xy[A, :]
|
||||
self._edge = np.r_[Utils.mkvc(length2D(edge1)), Utils.mkvc(length2D(edge2))]
|
||||
self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge]
|
||||
self._edge = np.r_[Utils.mkvc(length2D(edge1)),
|
||||
Utils.mkvc(length2D(edge2))]
|
||||
self._tangents = np.r_[edge1, edge2]/np.c_[self._edge,
|
||||
self._edge]
|
||||
elif(self.dim == 3):
|
||||
xyz = self.gridN
|
||||
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy, self.nNz]))
|
||||
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
|
||||
self.nNy,
|
||||
self.nNz]))
|
||||
edge1 = xyz[D, :] - xyz[A, :]
|
||||
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy, self.nNz]))
|
||||
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx,
|
||||
self.nCy,
|
||||
self.nNz]))
|
||||
edge2 = xyz[B, :] - xyz[A, :]
|
||||
A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx, self.nNy, self.nCz]))
|
||||
A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx,
|
||||
self.nNy,
|
||||
self.nCz]))
|
||||
edge3 = xyz[E, :] - xyz[A, :]
|
||||
self._edge = np.r_[Utils.mkvc(length3D(edge1)), Utils.mkvc(length3D(edge2)), Utils.mkvc(length3D(edge3))]
|
||||
self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge]
|
||||
self._edge = np.r_[Utils.mkvc(length3D(edge1)),
|
||||
Utils.mkvc(length3D(edge2)),
|
||||
Utils.mkvc(length3D(edge3))]
|
||||
self._tangents = (np.r_[edge1, edge2, edge3] /
|
||||
np.c_[self._edge, self._edge, self._edge])
|
||||
return self._edge
|
||||
return self._edge
|
||||
return locals()
|
||||
_edge = None
|
||||
edge = property(**edge())
|
||||
|
||||
def tangents():
|
||||
doc = "Edge tangents."
|
||||
|
||||
def fget(self):
|
||||
if(self._tangents is None):
|
||||
@property
|
||||
def tangents(self):
|
||||
"""
|
||||
Edge tangents
|
||||
"""
|
||||
if getattr(self, '_tangents', None) is None:
|
||||
self.edge # calling .edge will create the tangents
|
||||
return self._tangents
|
||||
return locals()
|
||||
_tangents = None
|
||||
tangents = property(**tangents())
|
||||
|
||||
|
||||
|
||||
|
||||
+239
-176
@@ -18,13 +18,15 @@ def checkBC(bc):
|
||||
|
||||
for bc_i in bc:
|
||||
assert type(bc_i) is str, "each bc must be a string"
|
||||
assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'"
|
||||
assert bc_i in ['dirichlet', 'neumann'], ("each bc must be either,"
|
||||
"'dirichlet' or 'neumann'")
|
||||
return bc
|
||||
|
||||
|
||||
def ddxCellGrad(n, bc):
|
||||
"""
|
||||
Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1
|
||||
Create 1D derivative operator from cell-centers to nodes this means we
|
||||
go from n to n+1
|
||||
|
||||
For Cell-Centered **Dirichlet**, use a ghost point::
|
||||
|
||||
@@ -52,7 +54,8 @@ def ddxCellGrad(n, bc):
|
||||
"""
|
||||
bc = checkBC(bc)
|
||||
|
||||
D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr")
|
||||
D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n,
|
||||
format="csr")
|
||||
# Set the first side
|
||||
if(bc[0] == 'dirichlet'):
|
||||
D[0, 0] = 2
|
||||
@@ -65,10 +68,11 @@ def ddxCellGrad(n, bc):
|
||||
D[-1, -1] = 0
|
||||
return D
|
||||
|
||||
|
||||
def ddxCellGradBC(n, bc):
|
||||
"""
|
||||
|
||||
Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1
|
||||
Create 1D derivative operator from cell-centers to nodes this means we
|
||||
go from n to n+1
|
||||
|
||||
For Cell-Centered **Dirichlet**, use a ghost point::
|
||||
|
||||
@@ -99,7 +103,7 @@ def ddxCellGradBC(n, bc):
|
||||
"""
|
||||
bc = checkBC(bc)
|
||||
|
||||
ij = (np.array([0, n]),np.array([0, 1]))
|
||||
ij = (np.array([0, n]), np.array([0, 1]))
|
||||
vals = np.zeros(2)
|
||||
|
||||
# Set the first side
|
||||
@@ -112,7 +116,7 @@ def ddxCellGradBC(n, bc):
|
||||
vals[1] = 2
|
||||
elif(bc[1] == 'neumann'):
|
||||
vals[1] = 0
|
||||
D = sp.csr_matrix((vals, ij), shape=(n+1,2))
|
||||
D = sp.csr_matrix((vals, ij), shape=(n+1, 2))
|
||||
return D
|
||||
|
||||
|
||||
@@ -121,14 +125,16 @@ class DiffOperators(object):
|
||||
Class creates the differential operators that you need!
|
||||
"""
|
||||
def __init__(self):
|
||||
raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.')
|
||||
raise Exception('DiffOperators is a base class providing differential'
|
||||
'operators on meshes and cannot run on its own.'
|
||||
'Inherit to your favorite Mesh class.')
|
||||
|
||||
def faceDiv():
|
||||
doc = "Construct divergence operator (face-stg to cell-centres)."
|
||||
|
||||
def fget(self):
|
||||
if(self._faceDiv is None):
|
||||
# The number of cell centers in each direction
|
||||
@property
|
||||
def faceDiv(self):
|
||||
"""
|
||||
Construct divergence operator (face-stg to cell-centres).
|
||||
"""
|
||||
if getattr(self, '_faceDiv', None) is None:
|
||||
n = self.vnC
|
||||
# Compute faceDivergence operator on faces
|
||||
if(self.dim == 1):
|
||||
@@ -146,17 +152,15 @@ class DiffOperators(object):
|
||||
S = self.area
|
||||
V = self.vol
|
||||
self._faceDiv = sdiag(1/V)*D*sdiag(S)
|
||||
|
||||
return self._faceDiv
|
||||
return locals()
|
||||
_faceDiv = None
|
||||
faceDiv = property(**faceDiv())
|
||||
|
||||
def faceDivx():
|
||||
doc = "Construct divergence operator in the x component (face-stg to cell-centres)."
|
||||
|
||||
def fget(self):
|
||||
if(self._faceDivx is None):
|
||||
@property
|
||||
def faceDivx(self):
|
||||
"""
|
||||
Construct divergence operator in the x component (face-stg to
|
||||
cell-centres).
|
||||
"""
|
||||
if getattr(self, '_faceDivx', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
# Compute faceDivergence operator on faces
|
||||
@@ -172,16 +176,12 @@ class DiffOperators(object):
|
||||
self._faceDivx = sdiag(1/V)*D1*sdiag(S)
|
||||
|
||||
return self._faceDivx
|
||||
return locals()
|
||||
_faceDivx = None
|
||||
faceDivx = property(**faceDivx())
|
||||
|
||||
def faceDivy():
|
||||
doc = "Construct divergence operator in the y component (face-stg to cell-centres)."
|
||||
|
||||
def fget(self):
|
||||
if(self.dim < 2): return None
|
||||
if(self._faceDivy is None):
|
||||
@property
|
||||
def faceDivy(self):
|
||||
if(self.dim < 2):
|
||||
return None
|
||||
if getattr(self, '_faceDivy', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
# Compute faceDivergence operator on faces
|
||||
@@ -193,18 +193,17 @@ class DiffOperators(object):
|
||||
S = self.r(self.area, 'F', 'Fy', 'V')
|
||||
V = self.vol
|
||||
self._faceDivy = sdiag(1/V)*D2*sdiag(S)
|
||||
|
||||
return self._faceDivy
|
||||
return locals()
|
||||
_faceDivy = None
|
||||
faceDivy = property(**faceDivy())
|
||||
|
||||
def faceDivz():
|
||||
doc = "Construct divergence operator in the z component (face-stg to cell-centres)."
|
||||
|
||||
def fget(self):
|
||||
if(self.dim < 3): return None
|
||||
if(self._faceDivz is None):
|
||||
@property
|
||||
def faceDivz(self):
|
||||
"""
|
||||
Construct divergence operator in the z component (face-stg to
|
||||
cell-centres).
|
||||
"""
|
||||
if(self.dim < 3):
|
||||
return None
|
||||
if getattr(self, '_faceDivz', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
# Compute faceDivergence operator on faces
|
||||
@@ -213,17 +212,14 @@ class DiffOperators(object):
|
||||
S = self.r(self.area, 'F', 'Fz', 'V')
|
||||
V = self.vol
|
||||
self._faceDivz = sdiag(1/V)*D3*sdiag(S)
|
||||
|
||||
return self._faceDivz
|
||||
return locals()
|
||||
_faceDivz = None
|
||||
faceDivz = property(**faceDivz())
|
||||
|
||||
def nodalGrad():
|
||||
doc = "Construct gradient operator (nodes to edges)."
|
||||
|
||||
def fget(self):
|
||||
if(self._nodalGrad is None):
|
||||
@property
|
||||
def nodalGrad(self):
|
||||
"""
|
||||
Construct gradient operator (nodes to edges).
|
||||
"""
|
||||
if getattr(self, '_nodalGrad', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
# Compute divergence operator on faces
|
||||
@@ -242,15 +238,13 @@ class DiffOperators(object):
|
||||
L = self.edge
|
||||
self._nodalGrad = sdiag(1/L)*G
|
||||
return self._nodalGrad
|
||||
return locals()
|
||||
_nodalGrad = None
|
||||
nodalGrad = property(**nodalGrad())
|
||||
|
||||
def nodalLaplacian():
|
||||
doc = "Construct laplacian operator (nodes to edges)."
|
||||
|
||||
def fget(self):
|
||||
if(self._nodalLaplacian is None):
|
||||
@property
|
||||
def nodalLaplacian(self):
|
||||
"""
|
||||
Construct laplacian operator (nodes to edges).
|
||||
"""
|
||||
if getattr(self, '_nodalLaplacian', None) is None:
|
||||
print 'Warning: Laplacian has not been tested rigorously.'
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
@@ -274,22 +268,23 @@ class DiffOperators(object):
|
||||
L = L1 + L2 + L3
|
||||
self._nodalLaplacian = L
|
||||
return self._nodalLaplacian
|
||||
return locals()
|
||||
_nodalLaplacian = None
|
||||
nodalLaplacian = property(**nodalLaplacian())
|
||||
|
||||
def setCellGradBC(self, BC):
|
||||
"""
|
||||
Function that sets the boundary conditions for cell-centred derivative operators.
|
||||
Function that sets the boundary conditions for cell-centred derivative
|
||||
operators.
|
||||
|
||||
Examples::
|
||||
# Neumann in all directions
|
||||
BC = 'neumann'
|
||||
|
||||
BC = 'neumann' # Neumann in all directions
|
||||
BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else
|
||||
BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
|
||||
# Dirichlet else
|
||||
# 3D, Dirichlet in y Neumann else
|
||||
BC = ['neumann', 'dirichlet', 'neumann']
|
||||
|
||||
# 3D, Neumann in x on bottom of domain, Dirichlet else
|
||||
BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet']
|
||||
"""
|
||||
|
||||
if(type(BC) is str):
|
||||
BC = [BC]*self.dim
|
||||
if(type(BC) is list):
|
||||
@@ -323,26 +318,24 @@ class DiffOperators(object):
|
||||
G = sp.vstack((G1, G2, G3), format="csr")
|
||||
return G
|
||||
|
||||
def cellGrad():
|
||||
doc = "The cell centered Gradient, takes you to cell faces."
|
||||
|
||||
def fget(self):
|
||||
if(self._cellGrad is None):
|
||||
@property
|
||||
def cellGrad(self):
|
||||
"""
|
||||
The cell centered Gradient, takes you to cell faces.
|
||||
"""
|
||||
if getattr(self, '_cellGrad', None) is None:
|
||||
G = self._cellGradStencil()
|
||||
# Compute areas of cell faces & volumes
|
||||
S = self.area
|
||||
S = self.area # Compute areas of cell faces & volumes
|
||||
V = self.aveCC2F*self.vol # Average volume between adjacent cells
|
||||
self._cellGrad = sdiag(S/V)*G
|
||||
return self._cellGrad
|
||||
return locals()
|
||||
_cellGrad = None
|
||||
cellGrad = property(**cellGrad())
|
||||
|
||||
def cellGradBC():
|
||||
doc = "The cell centered Gradient boundary condition matrix"
|
||||
|
||||
def fget(self):
|
||||
if(self._cellGradBC is None):
|
||||
@property
|
||||
def cellGradBC(self):
|
||||
"""
|
||||
The cell centered Gradient boundary condition matrix
|
||||
"""
|
||||
if getattr(self, '_cellGradBC', None) is None:
|
||||
BC = self.setCellGradBC(self._cellGradBC_list)
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
@@ -361,9 +354,33 @@ class DiffOperators(object):
|
||||
V = self.aveCC2F*self.vol # Average volume between adjacent cells
|
||||
self._cellGradBC = sdiag(S/V)*G
|
||||
return self._cellGradBC
|
||||
return locals()
|
||||
_cellGradBC = None
|
||||
cellGradBC = property(**cellGradBC())
|
||||
|
||||
# def cellGradBC():
|
||||
# doc = "The cell centered Gradient boundary condition matrix"
|
||||
|
||||
# def fget(self):
|
||||
# if(self._cellGradBC is None):
|
||||
# BC = self.setCellGradBC(self._cellGradBC_list)
|
||||
# n = self.vnC
|
||||
# if(self.dim == 1):
|
||||
# G = ddxCellGradBC(n[0], BC[0])
|
||||
# elif(self.dim == 2):
|
||||
# G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
|
||||
# G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
|
||||
# G = sp.block_diag((G1, G2), format="csr")
|
||||
# elif(self.dim == 3):
|
||||
# G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
|
||||
# G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
|
||||
# G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
|
||||
# G = sp.block_diag((G1, G2, G3), format="csr")
|
||||
# # Compute areas of cell faces & volumes
|
||||
# S = self.area
|
||||
# V = self.aveCC2F*self.vol # Average volume between adjacent cells
|
||||
# self._cellGradBC = sdiag(S/V)*G
|
||||
# return self._cellGradBC
|
||||
# return locals()
|
||||
# _cellGradBC = None
|
||||
# cellGradBC = property(**cellGradBC())
|
||||
|
||||
def _cellGradxStencil(self):
|
||||
BC = ['neumann', 'neumann']
|
||||
@@ -376,11 +393,12 @@ class DiffOperators(object):
|
||||
G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
|
||||
return G1
|
||||
|
||||
|
||||
def cellGradx():
|
||||
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
|
||||
|
||||
def fget(self):
|
||||
@property
|
||||
def cellGradx(self):
|
||||
"""
|
||||
Cell centered Gradient in the x dimension. Has neumann boundary
|
||||
conditions.
|
||||
"""
|
||||
if getattr(self, '_cellGradx', None) is None:
|
||||
G1 = self._cellGradxStencil()
|
||||
# Compute areas of cell faces & volumes
|
||||
@@ -388,8 +406,6 @@ class DiffOperators(object):
|
||||
L = self.r(self.area/V, 'F','Fx', 'V')
|
||||
self._cellGradx = sdiag(L)*G1
|
||||
return self._cellGradx
|
||||
return locals()
|
||||
cellGradx = property(**cellGradx())
|
||||
|
||||
def _cellGradyStencil(self):
|
||||
if self.dim < 2: return None
|
||||
@@ -401,19 +417,17 @@ class DiffOperators(object):
|
||||
G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
|
||||
return G2
|
||||
|
||||
def cellGrady():
|
||||
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
|
||||
def fget(self):
|
||||
if self.dim < 2: return None
|
||||
@property
|
||||
def cellGrady(self):
|
||||
if self.dim < 2:
|
||||
return None
|
||||
if getattr(self, '_cellGrady', None) is None:
|
||||
G2 = self._cellGradyStencil()
|
||||
# Compute areas of cell faces & volumes
|
||||
V = self.aveCC2F*self.vol
|
||||
L = self.r(self.area/V, 'F','Fy', 'V')
|
||||
L = self.r(self.area/V, 'F', 'Fy', 'V')
|
||||
self._cellGrady = sdiag(L)*G2
|
||||
return self._cellGrady
|
||||
return locals()
|
||||
cellGrady = property(**cellGrady())
|
||||
|
||||
def _cellGradzStencil(self):
|
||||
if self.dim < 3: return None
|
||||
@@ -422,34 +436,33 @@ class DiffOperators(object):
|
||||
G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
|
||||
return G3
|
||||
|
||||
def cellGradz():
|
||||
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
|
||||
def fget(self):
|
||||
if self.dim < 3: return None
|
||||
@property
|
||||
def cellGradz(self):
|
||||
"""
|
||||
Cell centered Gradient in the x dimension. Has neumann boundary
|
||||
conditions.
|
||||
"""
|
||||
if self.dim < 3:
|
||||
return None
|
||||
if getattr(self, '_cellGradz', None) is None:
|
||||
G3 = self._cellGradzStencil()
|
||||
# Compute areas of cell faces & volumes
|
||||
V = self.aveCC2F*self.vol
|
||||
L = self.r(self.area/V, 'F','Fz', 'V')
|
||||
L = self.r(self.area/V, 'F', 'Fz', 'V')
|
||||
self._cellGradz = sdiag(L)*G3
|
||||
return self._cellGradz
|
||||
return locals()
|
||||
cellGradz = property(**cellGradz())
|
||||
|
||||
def edgeCurl():
|
||||
doc = "Construct the 3D curl operator."
|
||||
|
||||
def fget(self):
|
||||
if(self._edgeCurl is None):
|
||||
@property
|
||||
def edgeCurl(self):
|
||||
"""
|
||||
Construct the 3D curl operator.
|
||||
"""
|
||||
if getattr(self, '_edgeCurl', None) is None:
|
||||
assert self.dim > 1, "Edge Curl only programed for 2 or 3D."
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
|
||||
# Compute lengths of cell edges
|
||||
L = self.edge
|
||||
|
||||
# Compute areas of cell faces
|
||||
S = self.area
|
||||
n = self.vnC # The number of cell centers in each direction
|
||||
L = self.edge # Compute lengths of cell edges
|
||||
S = self.area # Compute areas of cell faces
|
||||
|
||||
# Compute divergence operator on faces
|
||||
if self.dim == 2:
|
||||
@@ -477,11 +490,7 @@ class DiffOperators(object):
|
||||
sp.hstack((-D21, D12, O3))), format="csr")
|
||||
|
||||
self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
|
||||
|
||||
return self._edgeCurl
|
||||
return locals()
|
||||
_edgeCurl = None
|
||||
edgeCurl = property(**edgeCurl())
|
||||
|
||||
def getBCProjWF(self, BC, discretization='CC'):
|
||||
"""
|
||||
@@ -489,16 +498,19 @@ class DiffOperators(object):
|
||||
The weak form boundary condition projection matrices.
|
||||
|
||||
Examples::
|
||||
# Neumann in all directions
|
||||
BC = 'neumann'
|
||||
|
||||
BC = 'neumann' # Neumann in all directions
|
||||
BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else
|
||||
BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
|
||||
# Dirichlet else
|
||||
# 3D, Dirichlet in y Neumann else
|
||||
BC = ['neumann', 'dirichlet', 'neumann']
|
||||
|
||||
# 3D, Neumann in x on bottom of domain, Dirichlet else
|
||||
BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet']
|
||||
"""
|
||||
|
||||
if discretization is not 'CC':
|
||||
raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
|
||||
raise NotImplementedError('Boundary conditions only implemented'
|
||||
'for CC discretization.')
|
||||
|
||||
if(type(BC) is str):
|
||||
BC = [BC for _ in self.vnC] # Repeat the str self.dim times
|
||||
@@ -510,35 +522,34 @@ class DiffOperators(object):
|
||||
for i, bc_i in enumerate(BC):
|
||||
BC[i] = checkBC(bc_i)
|
||||
|
||||
|
||||
def projDirichlet(n, bc):
|
||||
bc = checkBC(bc)
|
||||
ij = ([0,n], [0,1])
|
||||
vals = [0,0]
|
||||
ij = ([0, n], [0, 1])
|
||||
vals = [0, 0]
|
||||
if(bc[0] == 'dirichlet'):
|
||||
vals[0] = -1
|
||||
if(bc[1] == 'dirichlet'):
|
||||
vals[1] = 1
|
||||
return sp.csr_matrix((vals, ij), shape=(n+1,2))
|
||||
return sp.csr_matrix((vals, ij), shape=(n+1, 2))
|
||||
|
||||
def projNeumannIn(n, bc):
|
||||
bc = checkBC(bc)
|
||||
P = sp.identity(n+1).tocsr()
|
||||
if(bc[0] == 'neumann'):
|
||||
P = P[1:,:]
|
||||
P = P[1:, :]
|
||||
if(bc[1] == 'neumann'):
|
||||
P = P[:-1,:]
|
||||
P = P[:-1, :]
|
||||
return P
|
||||
|
||||
def projNeumannOut(n, bc):
|
||||
bc = checkBC(bc)
|
||||
ij = ([0, 1],[0, n])
|
||||
ij = ([0, 1], [0, n])
|
||||
vals = [0,0]
|
||||
if(bc[0] == 'neumann'):
|
||||
vals[0] = 1
|
||||
if(bc[1] == 'neumann'):
|
||||
vals[1] = 1
|
||||
return sp.csr_matrix((vals, ij), shape=(2,n+1))
|
||||
return sp.csr_matrix((vals, ij), shape=(2, n+1))
|
||||
|
||||
n = self.vnC
|
||||
indF = self.faceBoundaryInd
|
||||
@@ -550,6 +561,7 @@ class DiffOperators(object):
|
||||
Pin = projNeumannIn(n[0], BC[0])
|
||||
|
||||
Pout = projNeumannOut(n[0], BC[0])
|
||||
|
||||
elif(self.dim == 2):
|
||||
Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
|
||||
Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
|
||||
@@ -564,12 +576,14 @@ class DiffOperators(object):
|
||||
P1 = sp.kron(speye(n[1]), projNeumannOut(n[0], BC[0]))
|
||||
P2 = sp.kron(projNeumannOut(n[1], BC[1]), speye(n[0]))
|
||||
Pout = sp.block_diag((P1, P2), format="csr")
|
||||
|
||||
elif(self.dim == 3):
|
||||
Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
|
||||
Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
|
||||
Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
|
||||
Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
|
||||
indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
|
||||
indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] |
|
||||
indF[5])]
|
||||
Pbc = Pbc*sdiag(self.area[indF])
|
||||
|
||||
P1 = kron3(speye(n[2]), speye(n[1]), projNeumannIn(n[0], BC[0]))
|
||||
@@ -586,36 +600,36 @@ class DiffOperators(object):
|
||||
|
||||
def getBCProjWF_simple(self, discretization='CC'):
|
||||
"""
|
||||
|
||||
The weak form boundary condition projection matrices
|
||||
when mixed boundary condition is used
|
||||
|
||||
|
||||
"""
|
||||
|
||||
if discretization is not 'CC':
|
||||
raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
|
||||
raise NotImplementedError('Boundary conditions only implemented'
|
||||
'for CC discretization.')
|
||||
|
||||
def projBC(n):
|
||||
ij = ([0,n], [0,1])
|
||||
vals = [0,0]
|
||||
ij = ([0, n], [0, 1])
|
||||
vals = [0, 0]
|
||||
vals[0] = 1
|
||||
vals[1] = 1
|
||||
return sp.csr_matrix((vals, ij), shape=(n+1,2))
|
||||
return sp.csr_matrix((vals, ij), shape=(n+1, 2))
|
||||
|
||||
def projDirichlet(n, bc):
|
||||
bc = checkBC(bc)
|
||||
ij = ([0,n], [0,1])
|
||||
vals = [0,0]
|
||||
ij = ([0, n], [0, 1])
|
||||
vals = [0, 0]
|
||||
if(bc[0] == 'dirichlet'):
|
||||
vals[0] = -1
|
||||
if(bc[1] == 'dirichlet'):
|
||||
vals[1] = 1
|
||||
return sp.csr_matrix((vals, ij), shape=(n+1,2))
|
||||
return sp.csr_matrix((vals, ij), shape=(n+1, 2))
|
||||
|
||||
BC = [['dirichlet','dirichlet'],['dirichlet','dirichlet'],['dirichlet','dirichlet']]
|
||||
BC = [['dirichlet', 'dirichlet'], ['dirichlet', 'dirichlet'],
|
||||
['dirichlet', 'dirichlet']]
|
||||
n = self.vnC
|
||||
indF = self.faceBoundaryInd
|
||||
|
||||
if(self.dim == 1):
|
||||
Pbc = projDirichlet(n[0], BC[0])
|
||||
B = projBC(n[0])
|
||||
@@ -653,9 +667,11 @@ class DiffOperators(object):
|
||||
if(self.dim == 1):
|
||||
return self.aveFx2CC
|
||||
elif(self.dim == 2):
|
||||
return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC), format="csr")
|
||||
return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC),
|
||||
format="csr")
|
||||
elif(self.dim == 3):
|
||||
return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
|
||||
return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC,
|
||||
self.aveFz2CC), format="csr")
|
||||
|
||||
@property
|
||||
def aveF2CCV(self):
|
||||
@@ -665,11 +681,16 @@ class DiffOperators(object):
|
||||
elif(self.dim == 2):
|
||||
return sp.block_diag((self.aveFx2CC, self.aveFy2CC), format="csr")
|
||||
elif(self.dim == 3):
|
||||
return sp.block_diag((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
|
||||
return sp.block_diag((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC),
|
||||
format="csr")
|
||||
|
||||
@property
|
||||
def aveFx2CC(self):
|
||||
"Construct the averaging operator on cell faces in the x direction to cell centers."
|
||||
"""
|
||||
Construct the averaging operator on cell faces in the x direction to
|
||||
cell centers.
|
||||
"""
|
||||
|
||||
if getattr(self, '_aveFx2CC', None) is None:
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
@@ -682,8 +703,12 @@ class DiffOperators(object):
|
||||
|
||||
@property
|
||||
def aveFy2CC(self):
|
||||
"Construct the averaging operator on cell faces in the y direction to cell centers."
|
||||
if self.dim < 2: return None
|
||||
"""
|
||||
Construct the averaging operator on cell faces in the y direction to
|
||||
cell centers.
|
||||
"""
|
||||
if self.dim < 2:
|
||||
return None
|
||||
if getattr(self, '_aveFy2CC', None) is None:
|
||||
n = self.vnC
|
||||
if(self.dim == 2):
|
||||
@@ -694,7 +719,10 @@ class DiffOperators(object):
|
||||
|
||||
@property
|
||||
def aveFz2CC(self):
|
||||
"Construct the averaging operator on cell faces in the z direction to cell centers."
|
||||
"""
|
||||
Construct the averaging operator on cell faces in the z direction to
|
||||
cell centers.
|
||||
"""
|
||||
if self.dim < 3: return None
|
||||
if getattr(self, '_aveFz2CC', None) is None:
|
||||
n = self.vnC
|
||||
@@ -711,12 +739,18 @@ class DiffOperators(object):
|
||||
if(self.dim == 1):
|
||||
self._aveCC2F = avExtrap(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])),
|
||||
sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr")
|
||||
self._aveCC2F = sp.vstack((sp.kron(speye(n[1]),
|
||||
avExtrap(n[0])),
|
||||
sp.kron(avExtrap(n[1]),
|
||||
speye(n[0]))), format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])),
|
||||
kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])),
|
||||
kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr")
|
||||
self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]),
|
||||
avExtrap(n[0])),
|
||||
kron3(speye(n[2]), avExtrap(n[1]),
|
||||
speye(n[0])),
|
||||
kron3(avExtrap(n[2]), speye(n[1]),
|
||||
speye(n[0]))),
|
||||
format="csr")
|
||||
return self._aveCC2F
|
||||
|
||||
@property
|
||||
@@ -727,7 +761,8 @@ class DiffOperators(object):
|
||||
elif(self.dim == 2):
|
||||
return 0.5*sp.hstack((self.aveEx2CC, self.aveEy2CC), format="csr")
|
||||
elif(self.dim == 3):
|
||||
return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
|
||||
return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC,
|
||||
self.aveEz2CC), format="csr")
|
||||
|
||||
@property
|
||||
def aveE2CCV(self):
|
||||
@@ -737,11 +772,15 @@ class DiffOperators(object):
|
||||
elif(self.dim == 2):
|
||||
return sp.block_diag((self.aveEx2CC, self.aveEy2CC), format="csr")
|
||||
elif(self.dim == 3):
|
||||
return sp.block_diag((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
|
||||
return sp.block_diag((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC),
|
||||
format="csr")
|
||||
|
||||
@property
|
||||
def aveEx2CC(self):
|
||||
"Construct the averaging operator on cell edges in the x direction to cell centers."
|
||||
"""
|
||||
Construct the averaging operator on cell edges in the x direction to
|
||||
cell centers.
|
||||
"""
|
||||
if getattr(self, '_aveEx2CC', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
@@ -755,8 +794,12 @@ class DiffOperators(object):
|
||||
|
||||
@property
|
||||
def aveEy2CC(self):
|
||||
"Construct the averaging operator on cell edges in the y direction to cell centers."
|
||||
if self.dim < 2: return None
|
||||
"""
|
||||
Construct the averaging operator on cell edges in the y direction to
|
||||
cell centers.
|
||||
"""
|
||||
if self.dim < 2:
|
||||
return None
|
||||
if getattr(self, '_aveEy2CC', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
@@ -768,8 +811,12 @@ class DiffOperators(object):
|
||||
|
||||
@property
|
||||
def aveEz2CC(self):
|
||||
"Construct the averaging operator on cell edges in the z direction to cell centers."
|
||||
if self.dim < 3: return None
|
||||
"""
|
||||
Construct the averaging operator on cell edges in the z direction to
|
||||
cell centers.
|
||||
"""
|
||||
if self.dim < 3:
|
||||
return None
|
||||
if getattr(self, '_aveEz2CC', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
@@ -793,7 +840,10 @@ class DiffOperators(object):
|
||||
|
||||
@property
|
||||
def aveN2E(self):
|
||||
"Construct the averaging operator on cell nodes to cell edges, keeping each dimension separate."
|
||||
"""
|
||||
Construct the averaging operator on cell nodes to cell edges, keeping
|
||||
each dimension separate.
|
||||
"""
|
||||
|
||||
if getattr(self, '_aveN2E', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
@@ -802,16 +852,24 @@ class DiffOperators(object):
|
||||
self._aveN2E = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])),
|
||||
sp.kron(av(n[1]), speye(n[0]+1))), format="csr")
|
||||
sp.kron(av(n[1]), speye(n[0]+1))),
|
||||
format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])),
|
||||
kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)),
|
||||
kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr")
|
||||
self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1),
|
||||
av(n[0])),
|
||||
kron3(speye(n[2]+1), av(n[1]),
|
||||
speye(n[0]+1)),
|
||||
kron3(av(n[2]), speye(n[1]+1),
|
||||
speye(n[0]+1))),
|
||||
format="csr")
|
||||
return self._aveN2E
|
||||
|
||||
@property
|
||||
def aveN2F(self):
|
||||
"Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate."
|
||||
"""
|
||||
Construct the averaging operator on cell nodes to cell faces, keeping
|
||||
each dimension separate.
|
||||
"""
|
||||
if getattr(self, '_aveN2F', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
@@ -819,9 +877,14 @@ class DiffOperators(object):
|
||||
self._aveN2F = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)),
|
||||
sp.kron(speye(n[1]+1), av(n[0]))), format="csr")
|
||||
sp.kron(speye(n[1]+1), av(n[0]))),
|
||||
format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)),
|
||||
kron3(av(n[2]), speye(n[1]+1), av(n[0])),
|
||||
kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr")
|
||||
self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]),
|
||||
speye(n[0]+1)),
|
||||
kron3(av(n[2]), speye(n[1]+1),
|
||||
av(n[0])),
|
||||
kron3(speye(n[2]+1), av(n[1]),
|
||||
av(n[0]))),
|
||||
format="csr")
|
||||
return self._aveN2F
|
||||
|
||||
@@ -421,7 +421,7 @@ class InnerProducts(object):
|
||||
def _getEdgePx(M):
|
||||
"""Returns a function for creating projection matrices"""
|
||||
def Px(xEdge):
|
||||
assert xEdge == 'eX0', 'xEdge = %s, not eX0' % xEdge
|
||||
assert xEdge == 'eX0', 'xEdge = {0!s}, not eX0'.format(xEdge)
|
||||
return sp.identity(M.nC)
|
||||
return Px
|
||||
|
||||
|
||||
@@ -198,11 +198,11 @@ class TensorMeshIO(object):
|
||||
"""
|
||||
assert mesh.dim == 3
|
||||
s = ''
|
||||
s += '%i %i %i\n' %tuple(mesh.vnC)
|
||||
s += '{0:d} {1:d} {2:d}\n'.format(*tuple(mesh.vnC))
|
||||
origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated.
|
||||
origin.dtype = float
|
||||
|
||||
s += '%.2f %.2f %.2f\n' %tuple(origin)
|
||||
s += '{0:.2f} {1:.2f} {2:.2f}\n'.format(*tuple(origin))
|
||||
s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
|
||||
s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
|
||||
s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1])
|
||||
|
||||
@@ -23,8 +23,8 @@ class BaseTensorMesh(BaseMesh):
|
||||
h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i)
|
||||
elif type(h_i) is list:
|
||||
h_i = Utils.meshTensor(h_i)
|
||||
assert isinstance(h_i, np.ndarray), ("h[%i] is not a numpy array." % i)
|
||||
assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i)
|
||||
assert isinstance(h_i, np.ndarray), ("h[{0:d}] is not a numpy array.".format(i))
|
||||
assert len(h_i.shape) == 1, ("h[{0:d}] must be a 1D numpy array.".format(i))
|
||||
h[i] = h_i[:] # make a copy.
|
||||
|
||||
x0 = np.zeros(len(h))
|
||||
@@ -41,7 +41,7 @@ class BaseTensorMesh(BaseMesh):
|
||||
elif x_i == 'N':
|
||||
x0[i] = -h_i.sum()
|
||||
else:
|
||||
raise Exception("x0[%i] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative." % i)
|
||||
raise Exception("x0[{0:d}] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative.".format(i))
|
||||
|
||||
if isinstance(self, BaseRectangularMesh):
|
||||
BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0)
|
||||
@@ -239,7 +239,7 @@ class BaseTensorMesh(BaseMesh):
|
||||
'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)
|
||||
raise Exception('Symmetric CylMesh does not support {0!s} interpolation, as this variable does not exist.'.format(locType))
|
||||
|
||||
loc = Utils.asArray_N_x_Dim(loc, self.dim)
|
||||
|
||||
|
||||
@@ -177,7 +177,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
|
||||
return l
|
||||
|
||||
def __str__(self):
|
||||
outStr = ' ---- %sTreeMesh ---- '%('Oc' if self.dim == 3 else 'Quad')
|
||||
outStr = ' ---- {0!s}TreeMesh ---- '.format(('Oc' if self.dim == 3 else 'Quad'))
|
||||
def printH(hx, outStr=''):
|
||||
i = -1
|
||||
while True:
|
||||
@@ -213,7 +213,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
|
||||
outStr += printH(self.hy, outStr='\n hy:')
|
||||
outStr += printH(self.hz, outStr='\n hz:')
|
||||
outStr += '\n nC: {0:d}'.format(self.nC)
|
||||
outStr += '\n Fill: %2.2f%%'%(self.fill*100)
|
||||
outStr += '\n Fill: {0:2.2f}%'.format((self.fill*100))
|
||||
return outStr
|
||||
|
||||
@property
|
||||
@@ -2210,7 +2210,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
|
||||
|
||||
ax.set_xlabel('y' if normal == 'X' else 'x')
|
||||
ax.set_ylabel('y' if normal == 'Z' else 'z')
|
||||
ax.set_title('Slice %d, %s = %4.2f' % (ind,normal,indLoc))
|
||||
ax.set_title('Slice {0:d}, {1!s} = {2:4.2f}'.format(ind, normal, indLoc))
|
||||
|
||||
if grid:
|
||||
_ = antiNormalInd
|
||||
@@ -2240,7 +2240,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
|
||||
if key < 0 : #Handle negative indices
|
||||
key += len( self )
|
||||
if key >= len( self ) :
|
||||
raise IndexError, "The index (%d) is out of range."%key
|
||||
raise IndexError, "The index ({0:d}) is out of range.".format(key)
|
||||
|
||||
self._numberCells() # no-op if numbered
|
||||
index = self._i2cc[key]
|
||||
|
||||
+8
-8
@@ -171,7 +171,7 @@ class TensorView(object):
|
||||
iz = ix + iy*nX
|
||||
if iz < self.nCz:
|
||||
ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad,
|
||||
'#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large')
|
||||
'#{0:.0f}'.format(iz),color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large')
|
||||
|
||||
ax.set_title(vType)
|
||||
if showIt: plt.show()
|
||||
@@ -221,10 +221,10 @@ class TensorView(object):
|
||||
vTypeOpts = ['CC', 'CCv','N','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez']
|
||||
|
||||
# Some user error checking
|
||||
assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
|
||||
assert vType in vTypeOpts, "vType must be in ['{0!s}']".format("','".join(vTypeOpts))
|
||||
assert self.dim == 3, 'Must be a 3D mesh. Use plotImage.'
|
||||
assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts)
|
||||
assert normal in normalOpts, "normal must be in ['%s']" % "','".join(normalOpts)
|
||||
assert view in viewOpts, "view must be in ['{0!s}']".format("','".join(viewOpts))
|
||||
assert normal in normalOpts, "normal must be in ['{0!s}']".format("','".join(normalOpts))
|
||||
assert type(grid) is bool, 'grid must be a boolean'
|
||||
|
||||
szSliceDim = getattr(self, 'nC'+normal.lower()) #: Size of the sliced dimension
|
||||
@@ -295,7 +295,7 @@ class TensorView(object):
|
||||
|
||||
ax.set_xlabel('y' if normal == 'X' else 'x')
|
||||
ax.set_ylabel('y' if normal == 'Z' else 'z')
|
||||
ax.set_title('Slice %d' % ind)
|
||||
ax.set_title('Slice {0:.0f}'.format(ind))
|
||||
return out
|
||||
|
||||
|
||||
@@ -316,11 +316,11 @@ class TensorView(object):
|
||||
vTypeOptsV = ['CCv','F','E']
|
||||
vTypeOpts = vTypeOptsCC + vTypeOptsV
|
||||
if view == 'vec':
|
||||
assert vType in vTypeOptsV, "vType must be in ['%s'] when view='vec'" % "','".join(vTypeOptsV)
|
||||
assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
|
||||
assert vType in vTypeOptsV, "vType must be in ['{0!s}'] when view='vec'".format("','".join(vTypeOptsV))
|
||||
assert vType in vTypeOpts, "vType must be in ['{0!s}']".format("','".join(vTypeOpts))
|
||||
|
||||
viewOpts = ['real','imag','abs','vec']
|
||||
assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts)
|
||||
assert view in viewOpts, "view must be in ['{0!s}']".format("','".join(viewOpts))
|
||||
|
||||
|
||||
if ax is None:
|
||||
|
||||
@@ -121,7 +121,7 @@ class Minimize(object):
|
||||
@callback.setter
|
||||
def callback(self, value):
|
||||
if self.callback is not None:
|
||||
print 'The callback on the %s Optimization was replaced.' % self.__name__
|
||||
print 'The callback on the {0!s} Optimization was replaced.'.format(self.__name__)
|
||||
self._callback = value
|
||||
|
||||
|
||||
@@ -855,7 +855,7 @@ class NewtonRoot(object):
|
||||
if self.comments and self.doLS: print '\tLinesearch:\n'
|
||||
# Enter Linesearch
|
||||
while True and self.doLS:
|
||||
if self.comments: print '\t\tResid: %e\n'%norm(rt)
|
||||
if self.comments: print '\t\tResid: {0:e}\n'.format(norm(rt))
|
||||
if norm(rt) <= norm(r) or norm(rt) < self.tol:
|
||||
break
|
||||
|
||||
@@ -873,7 +873,7 @@ class NewtonRoot(object):
|
||||
if norm(rt) < self.tol:
|
||||
break
|
||||
if self.iter > self.maxIter:
|
||||
print 'NewtonRoot stopped by maxIters (%d). norm: %4.4e' % (self.maxIter, norm(rt))
|
||||
print 'NewtonRoot stopped by maxIters ({0:d}). norm: {1:4.4e}'.format(self.maxIter, norm(rt))
|
||||
break
|
||||
|
||||
return x
|
||||
|
||||
+1
-1
@@ -49,7 +49,7 @@ class BaseProblem(object):
|
||||
|
||||
def pair(self, d):
|
||||
"""Bind a survey to this problem instance using pointers."""
|
||||
assert isinstance(d, self.surveyPair), "Data object must be an instance of a %s class."%(self.surveyPair.__name__)
|
||||
assert isinstance(d, self.surveyPair), "Data object must be an instance of a {0!s} class.".format((self.surveyPair.__name__))
|
||||
if d.ispaired:
|
||||
raise Exception("The survey object is already paired to a problem. Use survey.unpair()")
|
||||
self._survey = d
|
||||
|
||||
+28
-28
@@ -19,85 +19,85 @@ class Property(object):
|
||||
return getattr(self, '_propertyLink', None)
|
||||
@propertyLink.setter
|
||||
def propertyLink(self, value):
|
||||
assert type(value) is tuple and len(value) == 2 and type(value[0]) is str and issubclass(value[1], Maps.IdentityMap), 'Use format: ("%s", Maps.ReciprocalMap)'%self.name
|
||||
assert type(value) is tuple and len(value) == 2 and type(value[0]) is str and issubclass(value[1], Maps.IdentityMap), 'Use format: ("{0!s}", Maps.ReciprocalMap)'.format(self.name)
|
||||
self._propertyLink = value
|
||||
|
||||
def _getMapProperty(self):
|
||||
prop = self
|
||||
def fget(self):
|
||||
return getattr(self, '_%sMap'%prop.name, None)
|
||||
return getattr(self, '_{0!s}Map'.format(prop.name), None)
|
||||
def fset(self, val):
|
||||
if prop.propertyLink is not None:
|
||||
linkName, linkMap = prop.propertyLink
|
||||
assert getattr(self, '%sMap'%linkName, None) is None, 'Cannot set both sides of a linked property.'
|
||||
assert getattr(self, '{0!s}Map'.format(linkName), None) is None, 'Cannot set both sides of a linked property.'
|
||||
# TODO: Check if the mapping can be correct
|
||||
setattr(self, '_%sMap'%prop.name, val)
|
||||
setattr(self, '_{0!s}Map'.format(prop.name), val)
|
||||
return property(fget=fget, fset=fset, doc=prop.doc)
|
||||
|
||||
def _getIndexProperty(self):
|
||||
prop = self
|
||||
def fget(self):
|
||||
return getattr(self, '_%sIndex'%prop.name, slice(None))
|
||||
return getattr(self, '_{0!s}Index'.format(prop.name), slice(None))
|
||||
def fset(self, val):
|
||||
setattr(self, '_%sIndex'%prop.name, val)
|
||||
setattr(self, '_{0!s}Index'.format(prop.name), val)
|
||||
return property(fget=fget, fset=fset, doc=prop.doc)
|
||||
|
||||
def _getProperty(self):
|
||||
prop = self
|
||||
def fget(self):
|
||||
mapping = getattr(self, '%sMap'%prop.name)
|
||||
mapping = getattr(self, '{0!s}Map'.format(prop.name))
|
||||
if mapping is None and prop.propertyLink is None:
|
||||
return prop.defaultVal
|
||||
|
||||
if mapping is None and prop.propertyLink is not None:
|
||||
linkName, linkMapClass = prop.propertyLink
|
||||
linkMap = linkMapClass(None)
|
||||
if getattr(self, '%sMap'%linkName, None) is None:
|
||||
if getattr(self, '{0!s}Map'.format(linkName), None) is None:
|
||||
return prop.defaultVal
|
||||
m = getattr(self, '%s'%linkName)
|
||||
m = getattr(self, '{0!s}'.format(linkName))
|
||||
return linkMap * m
|
||||
|
||||
m = getattr(self, '%sModel'%prop.name)
|
||||
m = getattr(self, '{0!s}Model'.format(prop.name))
|
||||
return mapping * m
|
||||
return property(fget=fget)
|
||||
|
||||
def _getModelDerivProperty(self):
|
||||
prop = self
|
||||
def fget(self):
|
||||
mapping = getattr(self, '%sMap'%prop.name)
|
||||
mapping = getattr(self, '{0!s}Map'.format(prop.name))
|
||||
if mapping is None and prop.propertyLink is None:
|
||||
return None
|
||||
|
||||
if mapping is None and prop.propertyLink is not None:
|
||||
linkName, linkMapClass = prop.propertyLink
|
||||
linkedMap = getattr(self, '%sMap'%linkName)
|
||||
linkedMap = getattr(self, '{0!s}Map'.format(linkName))
|
||||
if linkedMap is None:
|
||||
return None
|
||||
linkMap = linkMapClass(None) * linkedMap
|
||||
m = getattr(self, '%sModel'%linkName)
|
||||
m = getattr(self, '{0!s}Model'.format(linkName))
|
||||
return linkMap.deriv( m )
|
||||
|
||||
m = getattr(self, '%sModel'%prop.name)
|
||||
m = getattr(self, '{0!s}Model'.format(prop.name))
|
||||
return mapping.deriv( m )
|
||||
return property(fget=fget)
|
||||
|
||||
def _getModelProperty(self):
|
||||
prop = self
|
||||
def fget(self):
|
||||
mapping = getattr(self, '%sMap'%prop.name)
|
||||
mapping = getattr(self, '{0!s}Map'.format(prop.name))
|
||||
if mapping is None:
|
||||
return None
|
||||
index = getattr(self.propMap, '%sIndex'%prop.name)
|
||||
index = getattr(self.propMap, '{0!s}Index'.format(prop.name))
|
||||
return self.vector[index]
|
||||
return property(fget=fget)
|
||||
|
||||
def _getModelProjProperty(self):
|
||||
prop = self
|
||||
def fget(self):
|
||||
mapping = getattr(self, '%sMap'%prop.name)
|
||||
mapping = getattr(self, '{0!s}Map'.format(prop.name))
|
||||
if mapping is None:
|
||||
return None
|
||||
inds = getattr(self.propMap, '%sIndex'%prop.name)
|
||||
inds = getattr(self.propMap, '{0!s}Index'.format(prop.name))
|
||||
if type(inds) is slice:
|
||||
inds = range(*inds.indices(self.nP))
|
||||
nI, nP = len(inds),self.nP
|
||||
@@ -107,7 +107,7 @@ class Property(object):
|
||||
def _getModelMapProperty(self):
|
||||
prop = self
|
||||
def fget(self):
|
||||
return getattr(self.propMap, '_%sMap'%prop.name, None)
|
||||
return getattr(self.propMap, '_{0!s}Map'.format(prop.name), None)
|
||||
return property(fget=fget)
|
||||
|
||||
|
||||
@@ -123,7 +123,7 @@ class PropModel(object):
|
||||
inds = []
|
||||
if getattr(self, '_nP', None) is None:
|
||||
for name in self.propMap._properties:
|
||||
index = getattr(self.propMap, '%sIndex'%name, None)
|
||||
index = getattr(self.propMap, '{0!s}Index'.format(name), None)
|
||||
if index is not None:
|
||||
if type(index) is slice:
|
||||
inds += range(*index.indices(len(self.vector)))
|
||||
@@ -163,9 +163,9 @@ class _PropMapMetaClass(type):
|
||||
if prop.defaultInvProp:
|
||||
defaultInvProps += [p]
|
||||
if prop.propertyLink is not None:
|
||||
assert prop.propertyLink[0] in _properties, "You can only link to things that exist: '%s' is trying to link to '%s'"%(prop.name, prop.propertyLink[0])
|
||||
assert prop.propertyLink[0] in _properties, "You can only link to things that exist: '{0!s}' is trying to link to '{1!s}'".format(prop.name, prop.propertyLink[0])
|
||||
if len(defaultInvProps) > 1:
|
||||
raise Exception('You have more than one default inversion property: %s' % defaultInvProps)
|
||||
raise Exception('You have more than one default inversion property: {0!s}'.format(defaultInvProps))
|
||||
|
||||
newClass = super(_PropMapMetaClass, cls).__new__(cls, name, bases, attrs)
|
||||
|
||||
@@ -223,7 +223,7 @@ class PropMap(object):
|
||||
type(m[0]) is str and
|
||||
m[0] in self._properties and
|
||||
isinstance(m[1], Maps.IdentityMap)
|
||||
for m in maps]), "Use signature: [%s]" % (', '.join(["('%s', %sMap)"%(p,p) for p in self._properties]))
|
||||
for m in maps]), "Use signature: [{0!s}]".format((', '.join(["('{0!s}', {1!s}Map)".format(p, p) for p in self._properties])))
|
||||
if slices is None:
|
||||
slices = dict()
|
||||
else:
|
||||
@@ -236,8 +236,8 @@ class PropMap(object):
|
||||
|
||||
nP = 0
|
||||
for name, mapping in maps:
|
||||
setattr(self, '%sMap'%name, mapping)
|
||||
setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP)))
|
||||
setattr(self, '{0!s}Map'.format(name), mapping)
|
||||
setattr(self, '{0!s}Index'.format(name), slices.get(name, slice(nP, nP + mapping.nP)))
|
||||
nP += mapping.nP
|
||||
self.nP = nP
|
||||
|
||||
@@ -250,12 +250,12 @@ class PropMap(object):
|
||||
|
||||
def clearMaps(self):
|
||||
for name in self._properties:
|
||||
setattr(self, '%sMap'%name, None)
|
||||
setattr(self, '%sIndex'%name, None)
|
||||
setattr(self, '{0!s}Map'.format(name), None)
|
||||
setattr(self, '{0!s}Index'.format(name), None)
|
||||
|
||||
def __call__(self, vec):
|
||||
return self.PropModel(self, vec)
|
||||
|
||||
def __contains__(self, val):
|
||||
activeMaps = [name for name in self._properties if getattr(self, '%sMap'%name) is not None]
|
||||
activeMaps = [name for name in self._properties if getattr(self, '{0!s}Map'.format(name)) is not None]
|
||||
return val in activeMaps
|
||||
|
||||
+6
-6
@@ -26,7 +26,7 @@ class BaseRx(object):
|
||||
def rxType(self, value):
|
||||
known = self.knownRxTypes
|
||||
if known is not None:
|
||||
assert value in known, "rxType must be in ['%s']" % ("', '".join(known))
|
||||
assert value in known, "rxType must be in ['{0!s}']".format(("', '".join(known)))
|
||||
self._rxType = value
|
||||
|
||||
@property
|
||||
@@ -125,7 +125,7 @@ class BaseSrc(object):
|
||||
def __init__(self, rxList, **kwargs):
|
||||
assert type(rxList) is list, 'rxList must be a list'
|
||||
for rx in rxList:
|
||||
assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxPair.__name__
|
||||
assert isinstance(rx, self.rxPair), 'rxList must be a {0!s}'.format(self.rxPair.__name__)
|
||||
assert len(set(rxList)) == len(rxList), 'The rxList must be unique'
|
||||
self.uid = str(uuid.uuid4())
|
||||
self.rxList = rxList
|
||||
@@ -227,7 +227,7 @@ class BaseSurvey(object):
|
||||
@srcList.setter
|
||||
def srcList(self, value):
|
||||
assert type(value) is list, 'srcList must be a list'
|
||||
assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of %s' % self.srcPair.__name__
|
||||
assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of {0!s}'.format(self.srcPair.__name__)
|
||||
assert len(set(value)) == len(value), 'The srcList must be unique'
|
||||
self._srcList = value
|
||||
self._sourceOrder = dict()
|
||||
@@ -238,10 +238,10 @@ class BaseSurvey(object):
|
||||
sources = [sources]
|
||||
for src in sources:
|
||||
if getattr(src,'uid',None) is None:
|
||||
raise KeyError('Source does not have a uid: %s'%str(src))
|
||||
raise KeyError('Source does not have a uid: {0!s}'.format(str(src)))
|
||||
inds = map(lambda src: self._sourceOrder.get(src.uid, None), sources)
|
||||
if None in inds:
|
||||
raise KeyError('Some of the sources specified are not in this survey. %s'%str(inds))
|
||||
raise KeyError('Some of the sources specified are not in this survey. {0!s}'.format(str(inds)))
|
||||
return inds
|
||||
|
||||
@property
|
||||
@@ -263,7 +263,7 @@ class BaseSurvey(object):
|
||||
def pair(self, p):
|
||||
"""Bind a problem to this survey instance using pointers"""
|
||||
assert hasattr(p, 'surveyPair'), "Problem must have an attribute 'surveyPair'."
|
||||
assert isinstance(self, p.surveyPair), "Problem requires survey object must be an instance of a %s class."%(p.surveyPair.__name__)
|
||||
assert isinstance(self, p.surveyPair), "Problem requires survey object must be an instance of a {0!s} class.".format((p.surveyPair.__name__))
|
||||
if p.ispaired:
|
||||
raise Exception("The problem object is already paired to a survey. Use prob.unpair()")
|
||||
self._prob = p
|
||||
|
||||
+8
-9
@@ -4,7 +4,6 @@ from SimPEG.Utils import mkvc, sdiag, diagEst
|
||||
from SimPEG import Utils
|
||||
from SimPEG.Mesh import TensorMesh, CurvilinearMesh, CylMesh
|
||||
from SimPEG.Mesh.TreeMesh import TreeMesh as Tree
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import unittest
|
||||
import inspect
|
||||
@@ -200,10 +199,10 @@ class OrderTest(unittest.TestCase):
|
||||
print '_____________________________________________'
|
||||
print ' h | error | e(i-1)/e(i) | order'
|
||||
print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~'
|
||||
print '%4i | %8.2e |' % (nc, err)
|
||||
print '{0:4d} | {1:8.2e} |'.format(nc, err)
|
||||
else:
|
||||
order.append(np.log(err/err_old)/np.log(max_h/max_h_old))
|
||||
print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1])
|
||||
print '{0:4d} | {1:8.2e} | {2:6.4f} | {3:6.4f}'.format(nc, err, err_old/err, order[-1])
|
||||
err_old = err
|
||||
max_h_old = max_h
|
||||
print '---------------------------------------------'
|
||||
@@ -258,8 +257,8 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
|
||||
Tests.checkDerivative(simplePass, np.random.randn(5))
|
||||
"""
|
||||
|
||||
print "%s checkDerivative %s" % ('='*20, '='*20)
|
||||
print "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n%s" % ('-'*57)
|
||||
print "{0!s} checkDerivative {1!s}".format('='*20, '='*20)
|
||||
print "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n{0!s}".format(('-'*57))
|
||||
|
||||
f0, J0 = fctn(x0)
|
||||
|
||||
@@ -290,7 +289,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
|
||||
|
||||
order0 = np.log10(E0[:-1]/E0[1:])
|
||||
order1 = np.log10(E1[:-1]/E1[1:])
|
||||
print " %d %1.2e %1.3e %1.3e %1.3f" % (i, h[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1])
|
||||
print " {0:d} {1:1.2e} {2:1.3e} {3:1.3e} {4:1.3f}".format(i, h[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1])
|
||||
|
||||
# Ensure we are about precision
|
||||
order0 = order0[E0[1:] > eps]
|
||||
@@ -302,10 +301,10 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
|
||||
passTest = belowTol or correctOrder
|
||||
|
||||
if passTest:
|
||||
print "%s PASS! %s" % ('='*25, '='*25)
|
||||
print "{0!s} PASS! {1!s}".format('='*25, '='*25)
|
||||
print happiness[np.random.randint(len(happiness))]+'\n'
|
||||
else:
|
||||
print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57)
|
||||
print "{0!s}\n{1!s} FAIL! {2!s}\n{3!s}".format('*'*57, '<'*25, '>'*25, '*'*57)
|
||||
print sadness[np.random.randint(len(sadness))]+'\n'
|
||||
|
||||
|
||||
@@ -314,7 +313,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
|
||||
ax = ax or plt.subplot(111)
|
||||
ax.loglog(h, E0, 'b')
|
||||
ax.loglog(h, E1, 'g--')
|
||||
ax.set_title('Check Derivative - %s' % ('PASSED :)' if passTest else 'FAILED :('))
|
||||
ax.set_title('Check Derivative - {0!s}'.format(('PASSED :)' if passTest else 'FAILED :(')))
|
||||
ax.set_xlabel('h')
|
||||
ax.set_ylabel('Error')
|
||||
leg = ax.legend(['$\mathcal{O}(h)$', '$\mathcal{O}(h^2)$'], loc='best',
|
||||
|
||||
@@ -8,7 +8,7 @@ def _checkAccuracy(A, b, X, accuracyTol):
|
||||
if nrm_b > 0:
|
||||
nrm /= nrm_b
|
||||
if nrm > accuracyTol:
|
||||
msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol)
|
||||
msg = '### SolverWarning ###: Accuracy on solve is above tolerance: {0:e} > {1:e}'.format(nrm, accuracyTol)
|
||||
print msg
|
||||
warnings.warn(msg, RuntimeWarning)
|
||||
|
||||
|
||||
+15
-15
@@ -32,7 +32,7 @@ def memProfileWrapper(towrap, *funNames):
|
||||
if hasattr(towrap,f):
|
||||
attrs[f] = profile(getattr(towrap,f))
|
||||
else:
|
||||
print '%s not found in %s Class' % (f, towrap.__name__)
|
||||
print '{0!s} not found in {1!s} Class'.format(f, towrap.__name__)
|
||||
|
||||
return type(towrap.__name__ + 'MemProfileWrap', (towrap,), attrs)
|
||||
|
||||
@@ -65,7 +65,7 @@ def setKwargs(obj, ignore=None, **kwargs):
|
||||
if hasattr(obj, attr):
|
||||
setattr(obj, attr, kwargs[attr])
|
||||
else:
|
||||
raise Exception('%s attr is not recognized' % attr)
|
||||
raise Exception('{0!s} attr is not recognized'.format(attr))
|
||||
|
||||
hook(obj,hook, silent=True)
|
||||
hook(obj,setKwargs, silent=True)
|
||||
@@ -74,7 +74,7 @@ def printTitles(obj, printers, name='Print Titles', pad=''):
|
||||
titles = ''
|
||||
widths = 0
|
||||
for printer in printers:
|
||||
titles += ('{:^%i}'%printer['width']).format(printer['title']) + ''
|
||||
titles += ('{{:^{0:d}}}'.format(printer['width'])).format(printer['title']) + ''
|
||||
widths += printer['width']
|
||||
print pad + "{0} {1} {0}".format('='*((widths-1-len(name))/2), name)
|
||||
print pad + titles
|
||||
@@ -83,7 +83,7 @@ def printTitles(obj, printers, name='Print Titles', pad=''):
|
||||
def printLine(obj, printers, pad=''):
|
||||
values = ''
|
||||
for printer in printers:
|
||||
values += ('{:^%i}'%printer['width']).format(printer['format'] % printer['value'](obj))
|
||||
values += ('{{:^{0:d}}}'.format(printer['width'])).format(printer['format'] % printer['value'](obj))
|
||||
print pad + values
|
||||
|
||||
def checkStoppers(obj, stoppers):
|
||||
@@ -104,12 +104,12 @@ def checkStoppers(obj, stoppers):
|
||||
return (len(optimal)>0 and all(optimal)) | (len(critical)>0 and any(critical))
|
||||
|
||||
def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'):
|
||||
print pad + "%s%s%s" % ('-'*25,stop,'-'*25)
|
||||
print pad + "{0!s}{1!s}{2!s}".format('-'*25, stop, '-'*25)
|
||||
for stopper in stoppers:
|
||||
l = stopper['left'](obj)
|
||||
r = stopper['right'](obj)
|
||||
print pad + stopper['str'] % (l<=r,l,r)
|
||||
print pad + "%s%s%s" % ('-'*25,done,'-'*25)
|
||||
print pad + "{0!s}{1!s}{2!s}".format('-'*25, done, '-'*25)
|
||||
|
||||
def callHooks(match, mainFirst=False):
|
||||
"""
|
||||
@@ -144,14 +144,14 @@ def callHooks(match, mainFirst=False):
|
||||
|
||||
|
||||
extra = """
|
||||
If you have things that also need to run in the method %s, you can create a method::
|
||||
If you have things that also need to run in the method {0!s}, you can create a method::
|
||||
|
||||
def _%s*(self, ... ):
|
||||
def _{1!s}*(self, ... ):
|
||||
pass
|
||||
|
||||
Where the * can be any string. If present, _%s* will be called at the start of the default %s call.
|
||||
Where the * can be any string. If present, _{2!s}* will be called at the start of the default {3!s} call.
|
||||
You may also completely overwrite this function.
|
||||
""" % (match, match, match, match)
|
||||
""".format(match, match, match, match)
|
||||
doc = wrapper.__doc__
|
||||
wrapper.__doc__ = ('' if doc is None else doc) + extra
|
||||
return wrapper
|
||||
@@ -186,7 +186,7 @@ def asArray_N_x_Dim(pts, dim):
|
||||
elif len(pts.shape) == 1:
|
||||
pts = pts[:,np.newaxis]
|
||||
|
||||
assert pts.shape[1] == dim, "pts must be a column vector of shape (nPts, %d) not (%d, %d)" % ((dim,)+pts.shape)
|
||||
assert pts.shape[1] == dim, "pts must be a column vector of shape (nPts, {0:d}) not ({1:d}, {2:d})".format(*((dim,)+pts.shape))
|
||||
|
||||
return pts
|
||||
|
||||
@@ -207,17 +207,17 @@ def requires(var):
|
||||
|
||||
.. note::
|
||||
|
||||
To use survey.%s(), SimPEG requires that a problem be bound to the survey.
|
||||
To use survey.{0!s}(), SimPEG requires that a problem be bound to the survey.
|
||||
If a problem has not been bound, an Exception will be raised.
|
||||
To bind a problem to the Data object::
|
||||
|
||||
survey.pair(myProblem)
|
||||
|
||||
""" % f.__name__
|
||||
""".format(f.__name__)
|
||||
else:
|
||||
extra = """
|
||||
To use *%s* method, SimPEG requires that the %s be specified.
|
||||
""" % (f.__name__, var)
|
||||
To use *{0!s}* method, SimPEG requires that the {1!s} be specified.
|
||||
""".format(f.__name__, var)
|
||||
@wraps(f)
|
||||
def requiresVarWrapper(self,*args,**kwargs):
|
||||
if getattr(self, var, None) is None:
|
||||
|
||||
@@ -80,7 +80,7 @@ def indexCube(nodes, gridSize, n=None):
|
||||
# Make sure that we choose from the possible nodes.
|
||||
possibleNodes = 'ABCD' if gridSize.size == 2 else 'ABCDEFGH'
|
||||
for node in nodes:
|
||||
assert node in possibleNodes, "Nodes must be chosen from: '%s'" % possibleNodes
|
||||
assert node in possibleNodes, "Nodes must be chosen from: '{0!s}'".format(possibleNodes)
|
||||
dim = gridSize.size
|
||||
if n is None:
|
||||
n = gridSize - 1
|
||||
|
||||
@@ -278,7 +278,7 @@ class TensorType(object):
|
||||
else:
|
||||
raise Exception('Unexpected shape of tensor')
|
||||
def __str__(self):
|
||||
return 'TensorType[%i]: %s' % (self._tt, self._tts)
|
||||
return 'TensorType[{0:d}]: {1!s}'.format(self._tt, self._tts)
|
||||
def __eq__(self, v): return self._tt == v
|
||||
def __le__(self, v): return self._tt <= v
|
||||
def __ge__(self, v): return self._tt >= v
|
||||
|
||||
@@ -26,7 +26,7 @@ def surface2ind_topo(mesh, topo, gridLoc='CC'):
|
||||
gridTopo = Ftopo(XY).reshape(mesh.vnN[:2], order='F')
|
||||
|
||||
if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']:
|
||||
raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType)
|
||||
raise NotImplementedError('Nodal surface2ind_topo not implemented for {0!s} mesh'.format(mesh._meshType))
|
||||
|
||||
Nz = mesh.vectorNz[1:] # TODO: this will only work for tensor meshes
|
||||
actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F')
|
||||
@@ -47,7 +47,7 @@ def surface2ind_topo(mesh, topo, gridLoc='CC'):
|
||||
|
||||
gridTopo = Ftopo(mesh.vectorNx)
|
||||
if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']:
|
||||
raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType)
|
||||
raise NotImplementedError('Nodal surface2ind_topo not implemented for {0!s} mesh'.format(mesh._meshType))
|
||||
|
||||
Ny = mesh.vectorNy[1:] # TODO: this will only work for tensor meshes
|
||||
actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F')
|
||||
|
||||
+1
-1
@@ -266,7 +266,7 @@ def _supress_nonlocal_image_warn(self, msg, node):
|
||||
from docutils.utils import get_source_line
|
||||
|
||||
if not msg.startswith('nonlocal image URI found:'):
|
||||
self._warnfunc(msg, '%s:%s' % get_source_line(node))
|
||||
self._warnfunc(msg, '{0!s}:{1!s}'.format(*get_source_line(node)))
|
||||
|
||||
supress_nonlocal_image_warn()
|
||||
|
||||
|
||||
@@ -66,7 +66,7 @@ Numpy and Matlab
|
||||
Lessons in Python
|
||||
-----------------
|
||||
|
||||
* `Software Carpentry <http://software-carpentry.org/v4/python/index.html>`_
|
||||
* `Software Carpentry <http://swcarpentry.github.io/python-novice-inflammation/>`_
|
||||
* `Introduction to NumPy and Matplotlib <http://www.youtube.com/watch?v=3Fp1zn5ao2M>`_
|
||||
|
||||
Editing Python
|
||||
|
||||
@@ -47,7 +47,16 @@ direct current (DC) resistivity and induced polarization (IP) geophysical proble
|
||||
DC resistivity survey
|
||||
=====================
|
||||
|
||||
Electrical resistivity of subsurface materials is measured by causing an electrical current to flow in the earth between one pair of electrodes while the voltage across a second pair of electrodes is measured. The result is an "apparent" resistivity which is a value representing the weighted average resistivity over a volume of the earth. Variations in this measurement are caused by variations in the soil, rock, and pore fluid electrical resistivity. Surveys require contact with the ground, so they can be labour intensive. Results are sometimes interpreted directly, but more commonly, 1D, 2D or 3D models are estimated using inversion procedures (`GPG <http://www.eos.ubc.ca/courses/eosc350/content/>`_).
|
||||
Electrical resistivity of subsurface materials is measured by causing an
|
||||
electrical current to flow in the earth between one pair of electrodes while
|
||||
the voltage across a second pair of electrodes is measured. The result is an
|
||||
"apparent" resistivity which is a value representing the weighted average
|
||||
resistivity over a volume of the earth. Variations in this measurement are
|
||||
caused by variations in the soil, rock, and pore fluid electrical resistivity.
|
||||
Surveys require contact with the ground, so they can be labour intensive.
|
||||
Results are sometimes interpreted directly, but more commonly, 1D, 2D or 3D
|
||||
models are estimated using inversion procedures (`GPG
|
||||
<http://gpg.geosci.xyz>`_).
|
||||
|
||||
|
||||
Background
|
||||
|
||||
+14
-6
@@ -16,16 +16,24 @@ SimPEG Documentation
|
||||
:target: https://github.com/simpeg/simpeg/blob/master/LICENSE
|
||||
:alt: BSD 3 clause license.
|
||||
|
||||
.. image:: https://img.shields.io/travis/simpeg/simpeg.svg
|
||||
:target: https://travis-ci.org/simpeg/simpeg?branch=master
|
||||
.. image:: https://api.travis-ci.org/simpeg/simpeg.svg?branch=master
|
||||
:target: https://travis-ci.org/simpeg/simpeg
|
||||
:alt: Travis CI build status
|
||||
|
||||
.. image:: https://img.shields.io/coveralls/simpeg/simpeg.svg
|
||||
:target: https://coveralls.io/r/simpeg/simpeg?branch=master
|
||||
:alt: Coverage status
|
||||
.. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
|
||||
:alt: gitter chat room at https://gitter.im/simpeg/simpeg
|
||||
:target: https://gitter.im/simpeg/simpeg
|
||||
|
||||
.. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg
|
||||
:target: https://codecov.io/gh/simpeg/simpeg
|
||||
:target: https://codecov.io/gh/simpeg/simpeg
|
||||
|
||||
.. image:: https://www.quantifiedcode.com/api/v1/project/933aa3decf444538aa432c8817169b6d/badge.svg
|
||||
:target: https://www.quantifiedcode.com/app/project/933aa3decf444538aa432c8817169b6d
|
||||
:alt: Code issues
|
||||
|
||||
.. image:: https://api.codacy.com/project/badge/Grade/4fc959a5294a418fa21fc7bc3b3aa078
|
||||
:target: https://www.codacy.com/app/lindseyheagy/simpeg?utm_source=github.com&utm_medium=referral&utm_content=simpeg/simpeg&utm_campaign=Badge_Grade
|
||||
:alt: codacy
|
||||
|
||||
|
||||
Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
|
||||
|
||||
+1
-1
@@ -54,7 +54,7 @@ class Images(webapp2.RequestHandler):
|
||||
class Redirect(webapp2.RequestHandler):
|
||||
def get(self):
|
||||
path = str(self.request.path).split(os.path.sep)[3:]
|
||||
self.redirect(('/%s'%os.path.sep.join(path)), permanent=True)
|
||||
self.redirect(('/{0!s}'.format(os.path.sep.join(path))), permanent=True)
|
||||
|
||||
|
||||
class MainPage(webapp2.RequestHandler):
|
||||
|
||||
+147
-108
@@ -1,14 +1,16 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import Mesh, Problem, Fields, Survey, Utils
|
||||
import numpy as np
|
||||
|
||||
|
||||
class FieldsTest(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
|
||||
x = np.linspace(5,10,3)
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0.])
|
||||
srcLoc = np.r_[0,0,0.]
|
||||
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
|
||||
[0, 0, -30])
|
||||
x = np.linspace(5, 10, 3)
|
||||
XYZ = Utils.ndgrid(x, x, np.r_[0.])
|
||||
srcLoc = np.r_[0., 0., 0.]
|
||||
rxList0 = Survey.BaseRx(XYZ, 'exi')
|
||||
Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
|
||||
rxList1 = Survey.BaseRx(XYZ, 'bxi')
|
||||
@@ -18,10 +20,13 @@ class FieldsTest(unittest.TestCase):
|
||||
rxList3 = Survey.BaseRx(XYZ, 'bxi')
|
||||
Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
|
||||
Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
|
||||
srcList = [Src0,Src1,Src2,Src3,Src4]
|
||||
srcList = [Src0, Src1, Src2, Src3, Src4]
|
||||
survey = Survey.BaseSurvey(srcList=srcList)
|
||||
self.D = Survey.Data(survey)
|
||||
self.F = Problem.Fields(mesh, survey, knownFields={'phi':'CC','e':'E','b':'F'}, dtype={"phi":float,"e":complex,"b":complex})
|
||||
self.F = Problem.Fields(mesh, survey, knownFields={'phi': 'CC',
|
||||
'e': 'E', 'b': 'F'},
|
||||
dtype={"phi": float, "e": complex,
|
||||
"b": complex})
|
||||
self.Src0 = Src0
|
||||
self.Src1 = Src1
|
||||
self.mesh = mesh
|
||||
@@ -38,83 +43,94 @@ class FieldsTest(unittest.TestCase):
|
||||
self.assertTrue('e' in F)
|
||||
|
||||
def test_overlappingFields(self):
|
||||
self.assertRaises(AssertionError, Problem.Fields, self.F.mesh, self.F.survey,
|
||||
knownFields={'b':'F'},
|
||||
aliasFields={'b':['b',(lambda F, b, ind: b)]})
|
||||
self.assertRaises(AssertionError, Problem.Fields, self.F.mesh,
|
||||
self.F.survey, knownFields={'b': 'F'},
|
||||
aliasFields={'b': ['b', (lambda F, b, ind: b)]})
|
||||
|
||||
def test_SetGet(self):
|
||||
F = self.F
|
||||
nSrc = F.survey.nSrc
|
||||
e = np.random.rand(F.mesh.nE, nSrc) + np.random.rand(F.mesh.nE, nSrc)*1j
|
||||
e = (np.random.rand(F.mesh.nE, nSrc) +
|
||||
np.random.rand(F.mesh.nE, nSrc)*1j)
|
||||
F[:, 'e'] = e
|
||||
b = np.random.rand(F.mesh.nF, nSrc) + np.random.rand(F.mesh.nF, nSrc)*1j
|
||||
b = (np.random.rand(F.mesh.nF, nSrc) +
|
||||
np.random.rand(F.mesh.nF, nSrc)*1j)
|
||||
F[:, 'b'] = b
|
||||
|
||||
self.assertTrue(np.all(F[:, 'e'] == e))
|
||||
self.assertTrue(np.all(F[:, 'b'] == b))
|
||||
F[:] = {'b':b,'e':e}
|
||||
F[:] = {'b': b, 'e': e}
|
||||
self.assertTrue(np.all(F[:, 'e'] == e))
|
||||
self.assertTrue(np.all(F[:, 'b'] == b))
|
||||
|
||||
for s in [0,0.0,np.r_[0],long(0)]:
|
||||
for s in [0, 0.0, np.r_[0], long(0)]:
|
||||
F[:, 'b'] = s
|
||||
self.assertTrue(np.all(F[:, 'b'] == b*0))
|
||||
|
||||
b = np.random.rand(F.mesh.nF,1)
|
||||
b = np.random.rand(F.mesh.nF, 1)
|
||||
F[self.Src0, 'b'] = b
|
||||
self.assertTrue(np.all(F[self.Src0, 'b'] == b))
|
||||
|
||||
b = np.random.rand(F.mesh.nF,1)
|
||||
b = np.random.rand(F.mesh.nF, 1)
|
||||
F[self.Src0, 'b'] = b
|
||||
self.assertTrue(np.all(F[self.Src0, 'b'] == b))
|
||||
|
||||
phi = np.random.rand(F.mesh.nC,2)
|
||||
F[[self.Src0,self.Src1], 'phi'] = phi
|
||||
self.assertTrue(np.all(F[[self.Src0,self.Src1], 'phi'] == phi))
|
||||
phi = np.random.rand(F.mesh.nC, 2)
|
||||
F[[self.Src0, self.Src1], 'phi'] = phi
|
||||
self.assertTrue(np.all(F[[self.Src0, self.Src1], 'phi'] == phi))
|
||||
|
||||
fdict = F[:,:]
|
||||
fdict = F[:, :]
|
||||
self.assertTrue(type(fdict) is dict)
|
||||
self.assertTrue(sorted([k for k in fdict]) == ['b','e','phi'])
|
||||
self.assertTrue(sorted([k for k in fdict]) == ['b', 'e', 'phi'])
|
||||
|
||||
b = np.random.rand(F.mesh.nF, 2)
|
||||
F[[self.Src0, self.Src1],'b'] = b
|
||||
self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,1))
|
||||
self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,1))
|
||||
self.assertTrue(np.all(F[self.Src0,'b'] == Utils.mkvc(b[:,0],2)))
|
||||
self.assertTrue(np.all(F[self.Src1,'b'] == Utils.mkvc(b[:,1],2)))
|
||||
F[[self.Src0, self.Src1], 'b'] = b
|
||||
self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF, 1))
|
||||
self.assertTrue(F[self.Src0, 'b'].shape == (F.mesh.nF, 1))
|
||||
self.assertTrue(np.all(F[self.Src0, 'b'] == Utils.mkvc(b[:, 0], 2)))
|
||||
self.assertTrue(np.all(F[self.Src1, 'b'] == Utils.mkvc(b[:, 1], 2)))
|
||||
|
||||
def test_assertions(self):
|
||||
freq = [self.Src0, self.Src1]
|
||||
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc)
|
||||
|
||||
def fun(): self.F[freq, 'b'] = bWrongSize
|
||||
self.assertRaises(ValueError, fun)
|
||||
|
||||
def fun(): self.F[-999.]
|
||||
self.assertRaises(KeyError, fun)
|
||||
|
||||
def fun(): self.F['notRight']
|
||||
self.assertRaises(KeyError, fun)
|
||||
def fun(): self.F[freq,'notThere']
|
||||
|
||||
def fun(): self.F[freq, 'notThere']
|
||||
self.assertRaises(KeyError, fun)
|
||||
|
||||
|
||||
class FieldsTest_Alias(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
|
||||
x = np.linspace(5,10,3)
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0.])
|
||||
srcLoc = np.r_[0,0,0.]
|
||||
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
|
||||
[0, 0, -30])
|
||||
x = np.linspace(5, 10, 3)
|
||||
XYZ = Utils.ndgrid(x, x, np.r_[0.])
|
||||
srcLoc = np.r_[0, 0, 0.]
|
||||
rxList0 = Survey.BaseRx(XYZ, 'exi')
|
||||
Src0 = Survey.BaseSrc([rxList0],loc=srcLoc)
|
||||
Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
|
||||
rxList1 = Survey.BaseRx(XYZ, 'bxi')
|
||||
Src1 = Survey.BaseSrc([rxList1],loc=srcLoc)
|
||||
Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
|
||||
rxList2 = Survey.BaseRx(XYZ, 'bxi')
|
||||
Src2 = Survey.BaseSrc([rxList2],loc=srcLoc)
|
||||
Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
|
||||
rxList3 = Survey.BaseRx(XYZ, 'bxi')
|
||||
Src3 = Survey.BaseSrc([rxList3],loc=srcLoc)
|
||||
Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3],loc=srcLoc)
|
||||
srcList = [Src0,Src1,Src2,Src3,Src4]
|
||||
Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
|
||||
Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
|
||||
srcList = [Src0, Src1, Src2, Src3, Src4]
|
||||
survey = Survey.BaseSurvey(srcList=srcList)
|
||||
self.F = Problem.Fields(mesh, survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',(lambda e, ind: self.F.mesh.edgeCurl * e)]})
|
||||
self.F = Problem.Fields(mesh, survey, knownFields={'e': 'E'},
|
||||
aliasFields={'b': ['e', 'F',
|
||||
(lambda e, ind:
|
||||
self.F.mesh.edgeCurl *
|
||||
e)]})
|
||||
self.Src0 = Src0
|
||||
self.Src1 = Src1
|
||||
self.mesh = mesh
|
||||
@@ -135,9 +151,9 @@ class FieldsTest_Alias(unittest.TestCase):
|
||||
nSrc = F.survey.nSrc
|
||||
e = np.random.rand(F.mesh.nE, nSrc)
|
||||
F[:, 'e'] = e
|
||||
self.assertTrue(np.all(F[:, 'b'] == F.mesh.edgeCurl * e ))
|
||||
self.assertTrue(np.all(F[:, 'b'] == F.mesh.edgeCurl * e))
|
||||
|
||||
e = np.random.rand(F.mesh.nE,1)
|
||||
e = np.random.rand(F.mesh.nE, 1)
|
||||
F[self.Src0, 'e'] = e
|
||||
self.assertTrue(np.all(F[self.Src0, 'b'] == F.mesh.edgeCurl * e))
|
||||
|
||||
@@ -149,19 +165,21 @@ class FieldsTest_Alias(unittest.TestCase):
|
||||
def alias(e, ind):
|
||||
self.assertTrue(ind[0] is self.Src0)
|
||||
return self.F.mesh.edgeCurl * e
|
||||
F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
|
||||
e = np.random.rand(F.mesh.nE,1)
|
||||
F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e': 'E'},
|
||||
aliasFields={'b': ['e', 'F', alias]})
|
||||
e = np.random.rand(F.mesh.nE, 1)
|
||||
F[self.Src0, 'e'] = e
|
||||
F[self.Src0, 'b']
|
||||
|
||||
|
||||
def alias(e, ind):
|
||||
self.assertTrue(type(ind) is list)
|
||||
self.assertTrue(ind[0] is self.Src0)
|
||||
self.assertTrue(ind[1] is self.Src1)
|
||||
return self.F.mesh.edgeCurl * e
|
||||
F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
|
||||
e = np.random.rand(F.mesh.nE,2)
|
||||
|
||||
F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e': 'E'},
|
||||
aliasFields={'b': ['e', 'F', alias]})
|
||||
e = np.random.rand(F.mesh.nE, 2)
|
||||
F[[self.Src0, self.Src1], 'e'] = e
|
||||
F[[self.Src0, self.Src1], 'b']
|
||||
|
||||
@@ -169,10 +187,11 @@ class FieldsTest_Alias(unittest.TestCase):
|
||||
class FieldsTest_Time(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
|
||||
x = np.linspace(5,10,3)
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0.])
|
||||
srcLoc = np.r_[0,0,0.]
|
||||
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
|
||||
[0, 0, -30])
|
||||
x = np.linspace(5, 10, 3)
|
||||
XYZ = Utils.ndgrid(x, x, np.r_[0.])
|
||||
srcLoc = np.r_[0, 0, 0.]
|
||||
rxList0 = Survey.BaseRx(XYZ, 'exi')
|
||||
Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
|
||||
rxList1 = Survey.BaseRx(XYZ, 'bxi')
|
||||
@@ -182,11 +201,13 @@ class FieldsTest_Time(unittest.TestCase):
|
||||
rxList3 = Survey.BaseRx(XYZ, 'bxi')
|
||||
Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
|
||||
Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
|
||||
srcList = [Src0,Src1,Src2,Src3,Src4]
|
||||
srcList = [Src0, Src1, Src2, Src3, Src4]
|
||||
survey = Survey.BaseSurvey(srcList=srcList)
|
||||
prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
|
||||
prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10., 3), (20., 2)])
|
||||
survey.pair(prob)
|
||||
self.F = Problem.TimeFields(mesh, survey, knownFields={'phi':'CC','e':'E','b':'F'})
|
||||
self.F = Problem.TimeFields(mesh, survey, knownFields={'phi': 'CC',
|
||||
'e': 'E',
|
||||
'b': 'F'})
|
||||
self.Src0 = Src0
|
||||
self.Src1 = Src1
|
||||
self.mesh = mesh
|
||||
@@ -216,81 +237,92 @@ class FieldsTest_Time(unittest.TestCase):
|
||||
|
||||
self.assertTrue(np.all(F[:, 'e'] == e))
|
||||
self.assertTrue(np.all(F[:, 'b'] == b))
|
||||
F[:] = {'b':b,'e':e}
|
||||
F[:] = {'b': b, 'e': e}
|
||||
self.assertTrue(np.all(F[:, 'e'] == e))
|
||||
self.assertTrue(np.all(F[:, 'b'] == b))
|
||||
|
||||
for s in [0,0.0,np.r_[0],long(0)]:
|
||||
for s in [0, 0.0, np.r_[0], long(0)]:
|
||||
F[:, 'b'] = s
|
||||
self.assertTrue(np.all(F[:, 'b'] == b*0))
|
||||
|
||||
b = np.random.rand(F.mesh.nF,1,nT)
|
||||
b = np.random.rand(F.mesh.nF, 1, nT)
|
||||
F[self.Src0, 'b'] = b
|
||||
self.assertTrue(np.all(F[self.Src0, 'b'] == b[:,0,:]))
|
||||
self.assertTrue(np.all(F[self.Src0, 'b'] == b[:, 0, :]))
|
||||
|
||||
b = np.random.rand(F.mesh.nF,1,nT)
|
||||
F[self.Src0, 'b', 0] = b[:,:,0]
|
||||
self.assertTrue(np.all(F[self.Src0, 'b', 0] == Utils.mkvc(b[:,0,0],2)))
|
||||
b = np.random.rand(F.mesh.nF, 1, nT)
|
||||
F[self.Src0, 'b', 0] = b[:, :, 0]
|
||||
self.assertTrue(np.all(F[self.Src0, 'b', 0] == Utils.mkvc(b[:, 0, 0],
|
||||
2)))
|
||||
|
||||
phi = np.random.rand(F.mesh.nC,2,nT)
|
||||
F[[self.Src0,self.Src1], 'phi'] = phi
|
||||
self.assertTrue(np.all(F[[self.Src0,self.Src1], 'phi'] == phi))
|
||||
phi = np.random.rand(F.mesh.nC, 2, nT)
|
||||
F[[self.Src0, self.Src1], 'phi'] = phi
|
||||
self.assertTrue(np.all(F[[self.Src0, self.Src1], 'phi'] == phi))
|
||||
|
||||
fdict = F[:]
|
||||
self.assertTrue(type(fdict) is dict)
|
||||
self.assertTrue(sorted([k for k in fdict]) == ['b','e','phi'])
|
||||
self.assertTrue(sorted([k for k in fdict]) == ['b', 'e', 'phi'])
|
||||
|
||||
b = np.random.rand(F.mesh.nF, 2, nT)
|
||||
F[[self.Src0, self.Src1],'b'] = b
|
||||
self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,nT))
|
||||
self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,nT))
|
||||
self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0,:]))
|
||||
self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1,:]))
|
||||
self.assertTrue(np.all(F[self.Src0,'b',1] == Utils.mkvc(b[:,0,1],2)))
|
||||
self.assertTrue(np.all(F[self.Src1,'b',1] == Utils.mkvc(b[:,1,1],2)))
|
||||
self.assertTrue(np.all(F[self.Src0,'b',4] == Utils.mkvc(b[:,0,4],2)))
|
||||
self.assertTrue(np.all(F[self.Src1,'b',4] == Utils.mkvc(b[:,1,4],2)))
|
||||
|
||||
F[[self.Src0, self.Src1], 'b'] = b
|
||||
self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF, nT))
|
||||
self.assertTrue(F[self.Src0, 'b'].shape == (F.mesh.nF, nT))
|
||||
self.assertTrue(np.all(F[self.Src0, 'b'] == b[:, 0, :]))
|
||||
self.assertTrue(np.all(F[self.Src1, 'b'] == b[:, 1, :]))
|
||||
self.assertTrue(np.all(F[self.Src0, 'b', 1] ==
|
||||
Utils.mkvc(b[:, 0, 1], 2)))
|
||||
self.assertTrue(np.all(F[self.Src1, 'b', 1] ==
|
||||
Utils.mkvc(b[:, 1, 1], 2)))
|
||||
self.assertTrue(np.all(F[self.Src0, 'b', 4] ==
|
||||
Utils.mkvc(b[:, 0, 4], 2)))
|
||||
self.assertTrue(np.all(F[self.Src1, 'b', 4] ==
|
||||
Utils.mkvc(b[:, 1, 4], 2)))
|
||||
|
||||
b = np.random.rand(F.mesh.nF, 2, nT)
|
||||
F[[self.Src0, self.Src1],'b', 0] = b[:,:,0]
|
||||
F[[self.Src0, self.Src1], 'b', 0] = b[:, :, 0]
|
||||
|
||||
def test_assertions(self):
|
||||
freq = [self.Src0, self.Src1]
|
||||
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc)
|
||||
|
||||
def fun(): self.F[freq, 'b'] = bWrongSize
|
||||
self.assertRaises(ValueError, fun)
|
||||
|
||||
def fun(): self.F[-999.]
|
||||
self.assertRaises(KeyError, fun)
|
||||
|
||||
def fun(): self.F['notRight']
|
||||
self.assertRaises(KeyError, fun)
|
||||
def fun(): self.F[freq,'notThere']
|
||||
|
||||
def fun(): self.F[freq, 'notThere']
|
||||
self.assertRaises(KeyError, fun)
|
||||
|
||||
|
||||
class FieldsTest_Time_Aliased(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
|
||||
x = np.linspace(5,10,3)
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0.])
|
||||
srcLoc = np.r_[0,0,0.]
|
||||
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
|
||||
[0, 0, -30])
|
||||
x = np.linspace(5, 10, 3)
|
||||
XYZ = Utils.ndgrid(x, x, np.r_[0.])
|
||||
srcLoc = np.r_[0, 0, 0.]
|
||||
rxList0 = Survey.BaseRx(XYZ, 'exi')
|
||||
Src0 = Survey.BaseSrc( [rxList0],loc=srcLoc)
|
||||
Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
|
||||
rxList1 = Survey.BaseRx(XYZ, 'bxi')
|
||||
Src1 = Survey.BaseSrc( [rxList1],loc=srcLoc)
|
||||
Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
|
||||
rxList2 = Survey.BaseRx(XYZ, 'bxi')
|
||||
Src2 = Survey.BaseSrc( [rxList2],loc=srcLoc)
|
||||
Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
|
||||
rxList3 = Survey.BaseRx(XYZ, 'bxi')
|
||||
Src3 = Survey.BaseSrc( [rxList3],loc=srcLoc)
|
||||
Src4 = Survey.BaseSrc( [rxList0, rxList1, rxList2, rxList3],loc=srcLoc)
|
||||
srcList = [Src0,Src1,Src2,Src3,Src4]
|
||||
Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
|
||||
Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
|
||||
srcList = [Src0, Src1, Src2, Src3, Src4]
|
||||
survey = Survey.BaseSurvey(srcList=srcList)
|
||||
prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
|
||||
prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10., 3), (20., 2)])
|
||||
survey.pair(prob)
|
||||
|
||||
def alias(b, srcInd, timeInd):
|
||||
return self.F.mesh.edgeCurl.T * b + timeInd
|
||||
self.F = Problem.TimeFields(mesh, survey, knownFields={'b':'F'}, aliasFields={'e':['b','E',alias]})
|
||||
self.F = Problem.TimeFields(mesh, survey, knownFields={'b': 'F'},
|
||||
aliasFields={'e': ['b', 'E', alias]})
|
||||
self.Src0 = Src0
|
||||
self.Src1 = Src1
|
||||
self.mesh = mesh
|
||||
@@ -307,32 +339,32 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
|
||||
self.assertTrue('e' in F)
|
||||
self.assertTrue('b' in F)
|
||||
|
||||
|
||||
def test_simpleAlias(self):
|
||||
F = self.F
|
||||
nSrc = F.survey.nSrc
|
||||
nT = F.survey.prob.nT + 1
|
||||
b = np.random.rand(F.mesh.nF, nSrc, nT)
|
||||
F[:, 'b', :] = b
|
||||
self.assertTrue(np.all(F[:, 'e', 0] == F.mesh.edgeCurl.T * b[:,:,0] ))
|
||||
self.assertTrue(np.all(F[:, 'e', 0] == F.mesh.edgeCurl.T * b[:, :, 0]))
|
||||
|
||||
e = range(nT)
|
||||
for i in range(nT):
|
||||
e[i] = F.mesh.edgeCurl.T*b[:,:,i] + i
|
||||
e[i] = e[i][:,:,np.newaxis]
|
||||
e[i] = F.mesh.edgeCurl.T*b[:, :, i] + i
|
||||
e[i] = e[i][:, :, np.newaxis]
|
||||
e = np.concatenate(e, axis=2)
|
||||
self.assertTrue(np.all(F[:, 'e', :] == e ))
|
||||
self.assertTrue(np.all(F[self.Src0, 'e', :] == e[:,0,:] ))
|
||||
self.assertTrue(np.all(F[self.Src1, 'e', :] == e[:,1,:] ))
|
||||
self.assertTrue(np.all(F[:, 'e', :] == e))
|
||||
self.assertTrue(np.all(F[self.Src0, 'e', :] == e[:, 0, :]))
|
||||
self.assertTrue(np.all(F[self.Src1, 'e', :] == e[:, 1, :]))
|
||||
for t in range(nT):
|
||||
self.assertTrue(np.all(F[self.Src1, 'e', t] == Utils.mkvc(e[:,1,t],2) ))
|
||||
self.assertTrue(np.all(F[self.Src1, 'e', t] ==
|
||||
Utils.mkvc(e[:, 1, t], 2)))
|
||||
|
||||
b = np.random.rand(F.mesh.nF,nT)
|
||||
F[self.Src0, 'b',:] = b
|
||||
b = np.random.rand(F.mesh.nF, nT)
|
||||
F[self.Src0, 'b', :] = b
|
||||
Cb = F.mesh.edgeCurl.T * b
|
||||
for i in range(Cb.shape[1]):
|
||||
Cb[:,i] += i
|
||||
self.assertTrue(np.all(F[self.Src0, 'e',:] == Cb))
|
||||
Cb[:, i] += i
|
||||
self.assertTrue(np.all(F[self.Src0, 'e', :] == Cb))
|
||||
|
||||
def f():
|
||||
F[self.Src0, 'e'] = F[self.Src0, 'e']
|
||||
@@ -341,35 +373,42 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
|
||||
def test_aliasFunction(self):
|
||||
nT = self.F.survey.prob.nT + 1
|
||||
count = [0]
|
||||
|
||||
def alias(e, srcInd, timeInd):
|
||||
count[0] += 1
|
||||
self.assertTrue(srcInd[0] is self.Src0)
|
||||
return self.F.mesh.edgeCurl * e
|
||||
F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
|
||||
e = np.random.rand(F.mesh.nE,1,nT)
|
||||
F = Problem.TimeFields(self.F.mesh, self.F.survey,
|
||||
knownFields={'e': 'E'},
|
||||
aliasFields={'b': ['e', 'F', alias]})
|
||||
e = np.random.rand(F.mesh.nE, 1, nT)
|
||||
F[self.Src0, 'e', :] = e
|
||||
F[self.Src0, 'b', :]
|
||||
self.assertTrue(count[0] == nT) # ensure that this is called for every time separately.
|
||||
e = np.random.rand(F.mesh.nE,1,1)
|
||||
# ensure that this is called for every time separately.
|
||||
self.assertTrue(count[0] == nT)
|
||||
e = np.random.rand(F.mesh.nE, 1, 1)
|
||||
F[self.Src0, 'e', 1] = e
|
||||
count[0] = 0
|
||||
F[self.Src0, 'b', 1]
|
||||
self.assertTrue(count[0] == 1) # ensure that this is called only once.
|
||||
|
||||
|
||||
def alias(e, srcInd, timeInd):
|
||||
count[0] += 1
|
||||
self.assertTrue(type(srcInd) is list)
|
||||
self.assertTrue(srcInd[0] is self.Src0)
|
||||
self.assertTrue(srcInd[1] is self.Src1)
|
||||
return self.F.mesh.edgeCurl * e
|
||||
F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
|
||||
e = np.random.rand(F.mesh.nE,2, nT)
|
||||
F = Problem.TimeFields(self.F.mesh, self.F.survey,
|
||||
knownFields={'e': 'E'},
|
||||
aliasFields={'b': ['e', 'F', alias]})
|
||||
e = np.random.rand(F.mesh.nE, 2, nT)
|
||||
F[[self.Src0, self.Src1], 'e', :] = e
|
||||
count[0] = 0
|
||||
F[[self.Src0, self.Src1], 'b', :]
|
||||
self.assertTrue(count[0] == nT) # ensure that this is called for every time separately.
|
||||
e = np.random.rand(F.mesh.nE,2, 1)
|
||||
|
||||
# ensure that this is called for every time separately.
|
||||
self.assertTrue(count[0] == nT)
|
||||
e = np.random.rand(F.mesh.nE, 2, 1)
|
||||
F[[self.Src0, self.Src1], 'e', 1] = e
|
||||
count[0] = 0
|
||||
F[[self.Src0, self.Src1], 'b', 1]
|
||||
|
||||
@@ -28,14 +28,14 @@ class RegularizationTests(unittest.TestCase):
|
||||
|
||||
for i, mesh in enumerate(self.meshlist):
|
||||
|
||||
print 'Testing %iD'%mesh.dim
|
||||
print 'Testing {0:d}D'.format(mesh.dim)
|
||||
|
||||
mapping = r.mapPair(mesh)
|
||||
reg = r(mesh, mapping=mapping)
|
||||
m = np.random.rand(mapping.nP)
|
||||
reg.mref = np.ones_like(m)*np.mean(m)
|
||||
|
||||
print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
|
||||
print 'Check: phi_m (mref) = {0:f}'.format(reg.eval(reg.mref))
|
||||
passed = reg.eval(reg.mref) < TOL
|
||||
self.assertTrue(passed)
|
||||
|
||||
@@ -56,7 +56,7 @@ class RegularizationTests(unittest.TestCase):
|
||||
|
||||
for i, mesh in enumerate(self.meshlist):
|
||||
|
||||
print 'Testing Active Cells %iD'%(mesh.dim)
|
||||
print 'Testing Active Cells {0:d}D'.format((mesh.dim))
|
||||
|
||||
if mesh.dim == 1:
|
||||
indActive = Utils.mkvc(mesh.gridCC <= 0.8)
|
||||
@@ -70,7 +70,7 @@ class RegularizationTests(unittest.TestCase):
|
||||
m = np.random.rand(mesh.nC)[indAct]
|
||||
reg.mref = np.ones_like(m)*np.mean(m)
|
||||
|
||||
print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
|
||||
print 'Check: phi_m (mref) = {0:f}'.format(reg.eval(reg.mref))
|
||||
passed = reg.eval(reg.mref) < TOL
|
||||
self.assertTrue(passed)
|
||||
|
||||
@@ -87,7 +87,7 @@ class RegularizationTests(unittest.TestCase):
|
||||
|
||||
for i, mesh in enumerate(self.meshlist):
|
||||
|
||||
print 'Testing %iD'%mesh.dim
|
||||
print 'Testing {0:d}D'.format(mesh.dim)
|
||||
|
||||
# mapping = r.mapPair(mesh)
|
||||
# reg = r(mesh, mapping=mapping)
|
||||
|
||||
+11
-11
@@ -14,9 +14,9 @@ class Doc_Test(unittest.TestCase):
|
||||
html_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build']+['html'])
|
||||
|
||||
check = subprocess.call(["sphinx-build", "-nW", "-b", "html", "-d",
|
||||
"%s"%(doctrees_path) ,
|
||||
"%s"%(self.path_to_docs),
|
||||
"%s"%(html_path)])
|
||||
"{0!s}".format((doctrees_path)) ,
|
||||
"{0!s}".format((self.path_to_docs)),
|
||||
"{0!s}".format((html_path))])
|
||||
assert check == 0
|
||||
|
||||
# def test_latex(self):
|
||||
@@ -29,15 +29,15 @@ class Doc_Test(unittest.TestCase):
|
||||
# "%s"%(latex_path)])
|
||||
# assert check == 0
|
||||
|
||||
# def test_linkcheck(self):
|
||||
# doctrees_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build']+['doctrees'])
|
||||
# link_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build'])
|
||||
def test_linkcheck(self):
|
||||
doctrees_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build']+['doctrees'])
|
||||
link_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build'])
|
||||
|
||||
# check = subprocess.call(["sphinx-build", "-nW", "-b", "linkcheck", "-d",
|
||||
# "%s"%(doctrees_path),
|
||||
# "%s"%(self.path_to_docs),
|
||||
# "%s"%(link_path)])
|
||||
# assert check == 0
|
||||
check = subprocess.call(["sphinx-build", "-nW", "-b", "linkcheck", "-d",
|
||||
"%s"%(doctrees_path),
|
||||
"%s"%(self.path_to_docs),
|
||||
"%s"%(link_path)])
|
||||
assert check == 0
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -21,7 +21,7 @@ SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVe
|
||||
|
||||
def adjointTest(fdemType, comp):
|
||||
prb = getFDEMProblem(fdemType, comp, SrcList, freq)
|
||||
print 'Adjoint %s formulation - %s' % (fdemType, comp)
|
||||
print 'Adjoint {0!s} formulation - {1!s}'.format(fdemType, comp)
|
||||
|
||||
m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
mu = np.ones(prb.mesh.nC)*MU
|
||||
|
||||
@@ -21,7 +21,7 @@ SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVe
|
||||
|
||||
def adjointTest(fdemType, comp):
|
||||
prb = getFDEMProblem(fdemType, comp, SrcList, freq)
|
||||
print 'Adjoint %s formulation - %s' % (fdemType, comp)
|
||||
print 'Adjoint {0!s} formulation - {1!s}'.format(fdemType, comp)
|
||||
|
||||
m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
mu = np.ones(prb.mesh.nC)*MU
|
||||
|
||||
@@ -26,7 +26,7 @@ SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVe
|
||||
def derivTest(fdemType, comp):
|
||||
|
||||
prb = getFDEMProblem(fdemType, comp, SrcType, freq)
|
||||
print '%s formulation - %s' % (fdemType, comp)
|
||||
print '{0!s} formulation - {1!s}'.format(fdemType, comp)
|
||||
x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
mu = np.log(np.ones(prb.mesh.nC)*MU)
|
||||
|
||||
|
||||
@@ -55,7 +55,7 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
|
||||
if showIt == True:
|
||||
plt.loglog(rx.times[bz_calc>0], bz_calc[bz_calc>0], 'r', rx.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--')
|
||||
plt.loglog(rx.times, abs(bz_ana), 'b*')
|
||||
plt.title('sig_half = %e'%sig_half)
|
||||
plt.title('sig_half = {0:e}'.format(sig_half))
|
||||
plt.show()
|
||||
|
||||
return log10diff
|
||||
|
||||
@@ -25,7 +25,7 @@ class compareInitFiles(unittest.TestCase):
|
||||
|
||||
def get(test):
|
||||
def test_func(self):
|
||||
print '\nTesting %s.run(plotIt=False)\n'%test
|
||||
print '\nTesting {0!s}.run(plotIt=False)\n'.format(test)
|
||||
getattr(Examples, test).run(plotIt=False)
|
||||
self.assertTrue(True)
|
||||
return test_func
|
||||
|
||||
@@ -121,7 +121,7 @@ class RichardsTests1D(unittest.TestCase):
|
||||
tol = TOL*(10**int(np.log10(np.abs(zJv))))
|
||||
passed = np.abs(vJz - zJv) < tol
|
||||
print 'Richards Adjoint Test - PressureHead'
|
||||
print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol)
|
||||
print '{0:4.4e} === {1:4.4e}, diff={2:4.4e} < {3:4e}'.format(vJz, zJv, np.abs(vJz - zJv), tol)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Sensitivity(self):
|
||||
@@ -193,7 +193,7 @@ class RichardsTests2D(unittest.TestCase):
|
||||
tol = TOL*(10**int(np.log10(np.abs(zJv))))
|
||||
passed = np.abs(vJz - zJv) < tol
|
||||
print '2D: Richards Adjoint Test - PressureHead'
|
||||
print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol)
|
||||
print '{0:4.4e} === {1:4.4e}, diff={2:4.4e} < {3:4e}'.format(vJz, zJv, np.abs(vJz - zJv), tol)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Sensitivity(self):
|
||||
@@ -265,7 +265,7 @@ class RichardsTests3D(unittest.TestCase):
|
||||
tol = TOL*(10**int(np.log10(np.abs(zJv))))
|
||||
passed = np.abs(vJz - zJv) < tol
|
||||
print '3D: Richards Adjoint Test - PressureHead'
|
||||
print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol)
|
||||
print '{0:4.4e} === {1:4.4e}, diff={2:4.4e} < {3:4e}'.format(vJz, zJv, np.abs(vJz - zJv), tol)
|
||||
self.assertTrue(passed,True)
|
||||
|
||||
def test_Sensitivity(self):
|
||||
|
||||
@@ -10,7 +10,9 @@ class BasicCurvTests(unittest.TestCase):
|
||||
a = np.array([1, 1, 1])
|
||||
b = np.array([1, 2])
|
||||
c = np.array([1, 4])
|
||||
gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h]
|
||||
|
||||
def gridIt(h): return [np.cumsum(np.r_[0, x]) for x in h]
|
||||
|
||||
X, Y = ndgrid(gridIt([a, b]), vector=False)
|
||||
self.TM2 = TensorMesh([a, b])
|
||||
self.Curv2 = CurvilinearMesh([X, Y])
|
||||
@@ -19,7 +21,10 @@ class BasicCurvTests(unittest.TestCase):
|
||||
self.Curv3 = CurvilinearMesh([X, Y, Z])
|
||||
|
||||
def test_area_3D(self):
|
||||
test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2])
|
||||
test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8,
|
||||
1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4,
|
||||
4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1,
|
||||
1, 2, 2, 2])
|
||||
self.assertTrue(np.all(self.Curv3.area == test_area))
|
||||
|
||||
def test_vol_3D(self):
|
||||
@@ -33,54 +38,85 @@ class BasicCurvTests(unittest.TestCase):
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_edge_3D(self):
|
||||
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
|
||||
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
|
||||
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2,
|
||||
2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2,
|
||||
2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4,
|
||||
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
|
||||
t1 = np.all(self.Curv3.edge == test_edge)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_edge_2D(self):
|
||||
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2])
|
||||
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2,
|
||||
2])
|
||||
t1 = np.all(self.Curv2.edge == test_edge)
|
||||
self.assertTrue(t1)
|
||||
|
||||
def test_tangents(self):
|
||||
T = self.Curv2.tangents
|
||||
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv2.nEx)))
|
||||
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv2.nEx)))
|
||||
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv2.nEy)))
|
||||
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv2.nEy)))
|
||||
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] ==
|
||||
np.ones(self.Curv2.nEx)))
|
||||
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] ==
|
||||
np.zeros(self.Curv2.nEx)))
|
||||
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] ==
|
||||
np.zeros(self.Curv2.nEy)))
|
||||
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] ==
|
||||
np.ones(self.Curv2.nEy)))
|
||||
|
||||
T = self.Curv3.tangents
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv3.nEx)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv3.nEx)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.Curv3.nEx)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] ==
|
||||
np.ones(self.Curv3.nEx)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] ==
|
||||
np.zeros(self.Curv3.nEx)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] ==
|
||||
np.zeros(self.Curv3.nEx)))
|
||||
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv3.nEy)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv3.nEy)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.Curv3.nEy)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] ==
|
||||
np.zeros(self.Curv3.nEy)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] ==
|
||||
np.ones(self.Curv3.nEy)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] ==
|
||||
np.zeros(self.Curv3.nEy)))
|
||||
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.Curv3.nEz)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.Curv3.nEz)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.Curv3.nEz)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] ==
|
||||
np.zeros(self.Curv3.nEz)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] ==
|
||||
np.zeros(self.Curv3.nEz)))
|
||||
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] ==
|
||||
np.ones(self.Curv3.nEz)))
|
||||
|
||||
def test_normals(self):
|
||||
N = self.Curv2.normals
|
||||
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv2.nFx)))
|
||||
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv2.nFx)))
|
||||
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv2.nFy)))
|
||||
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv2.nFy)))
|
||||
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] ==
|
||||
np.ones(self.Curv2.nFx)))
|
||||
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] ==
|
||||
np.zeros(self.Curv2.nFx)))
|
||||
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] ==
|
||||
np.zeros(self.Curv2.nFy)))
|
||||
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] ==
|
||||
np.ones(self.Curv2.nFy)))
|
||||
|
||||
N = self.Curv3.normals
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv3.nFx)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv3.nFx)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.Curv3.nFx)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] ==
|
||||
np.ones(self.Curv3.nFx)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] ==
|
||||
np.zeros(self.Curv3.nFx)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] ==
|
||||
np.zeros(self.Curv3.nFx)))
|
||||
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv3.nFy)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv3.nFy)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.Curv3.nFy)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] ==
|
||||
np.zeros(self.Curv3.nFy)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] ==
|
||||
np.ones(self.Curv3.nFy)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] ==
|
||||
np.zeros(self.Curv3.nFy)))
|
||||
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.Curv3.nFz)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.Curv3.nFz)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.Curv3.nFz)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] ==
|
||||
np.zeros(self.Curv3.nFz)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] ==
|
||||
np.zeros(self.Curv3.nFz)))
|
||||
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] ==
|
||||
np.ones(self.Curv3.nFz)))
|
||||
|
||||
def test_grid(self):
|
||||
self.assertTrue(np.all(self.Curv2.gridCC == self.TM2.gridCC))
|
||||
|
||||
@@ -2,12 +2,12 @@ import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import unittest
|
||||
import matplotlib.pyplot as plt
|
||||
from SimPEG import *
|
||||
from SimPEG import Mesh, Tests, Utils, Solver
|
||||
|
||||
MESHTYPES = ['uniformTensorMesh']
|
||||
|
||||
|
||||
def getxBCyBC_CC(mesh, alpha, beta, gamma):
|
||||
# def getxBCyBC(mesh, alpha, beta, gamma):
|
||||
"""
|
||||
This is a subfunction generating mixed-boundary condition:
|
||||
|
||||
@@ -17,17 +17,19 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
|
||||
|
||||
\rho \vec{j} = -\nabla \phi \phi
|
||||
|
||||
\alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega
|
||||
\alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r
|
||||
= \partial \Omega
|
||||
|
||||
xBC = f_1(\alpha, \beta, \gamma)
|
||||
yBC = f(\alpha, \beta, \gamma)
|
||||
|
||||
Computes xBC and yBC for cell-centered discretizations
|
||||
"""
|
||||
if mesh.dim == 1: #1D
|
||||
|
||||
if mesh.dim == 1: # 1D
|
||||
if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2):
|
||||
raise Exception("Lenght of list, alpha should be 2")
|
||||
fCCxm,fCCxp = mesh.cellBoundaryInd
|
||||
fCCxm, fCCxp = mesh.cellBoundaryInd
|
||||
nBC = fCCxm.sum()+fCCxp.sum()
|
||||
h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
|
||||
|
||||
@@ -50,11 +52,11 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
|
||||
xBC = np.r_[xBC_xm, xBC_xp]
|
||||
yBC = np.r_[yBC_xm, yBC_xp]
|
||||
|
||||
elif mesh.dim == 2: #2D
|
||||
elif mesh.dim == 2: # 2D
|
||||
if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4):
|
||||
raise Exception("Lenght of list, alpha should be 4")
|
||||
|
||||
fxm,fxp,fym,fyp = mesh.faceBoundaryInd
|
||||
fxm, fxp, fym, fyp = mesh.faceBoundaryInd
|
||||
nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
|
||||
|
||||
alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
|
||||
@@ -65,8 +67,10 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
|
||||
# h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
|
||||
# h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
|
||||
|
||||
h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
|
||||
h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
|
||||
h_xm = mesh.hx[0]*np.ones_like(alpha_xm)
|
||||
h_xp = mesh.hx[-1]*np.ones_like(alpha_xp)
|
||||
h_ym = mesh.hy[0]*np.ones_like(alpha_ym)
|
||||
h_yp = mesh.hy[-1]*np.ones_like(alpha_yp)
|
||||
|
||||
a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
|
||||
b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
|
||||
@@ -87,8 +91,10 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
|
||||
yBC_ym = 0.5*(1.-b_ym)
|
||||
yBC_yp = 0.5*(1.-1./b_yp)
|
||||
|
||||
sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
|
||||
sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
|
||||
sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm],
|
||||
np.arange(mesh.nFx)[fxp]])
|
||||
sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym],
|
||||
np.arange(mesh.nFy)[fyp]])
|
||||
|
||||
xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
|
||||
xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
|
||||
@@ -98,11 +104,11 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
|
||||
xBC = np.r_[xBC_x, xBC_y]
|
||||
yBC = np.r_[yBC_x, yBC_y]
|
||||
|
||||
elif mesh.dim == 3: #3D
|
||||
elif mesh.dim == 3: # 3D
|
||||
if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6):
|
||||
raise Exception("Lenght of list, alpha should be 6")
|
||||
# fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd
|
||||
fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd
|
||||
fxm, fxp, fym, fyp, fzm, fzp = mesh.faceBoundaryInd
|
||||
nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum()
|
||||
|
||||
alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
|
||||
@@ -116,9 +122,12 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
|
||||
# h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
|
||||
# h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2]
|
||||
|
||||
h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp)
|
||||
h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp)
|
||||
h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp)
|
||||
h_xm = mesh.hx[0]*np.ones_like(alpha_xm)
|
||||
h_xp = mesh.hx[-1]*np.ones_like(alpha_xp)
|
||||
h_ym = mesh.hy[0]*np.ones_like(alpha_ym)
|
||||
h_yp = mesh.hy[-1]*np.ones_like(alpha_yp)
|
||||
h_zm = mesh.hz[0]*np.ones_like(alpha_zm)
|
||||
h_zp = mesh.hz[-1]*np.ones_like(alpha_zp)
|
||||
|
||||
a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
|
||||
b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
|
||||
@@ -148,9 +157,12 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
|
||||
yBC_zm = 0.5*(1.-b_zm)
|
||||
yBC_zp = 0.5*(1.-1./b_zp)
|
||||
|
||||
sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
|
||||
sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
|
||||
sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]])
|
||||
sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm],
|
||||
np.arange(mesh.nFx)[fxp]])
|
||||
sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym],
|
||||
np.arange(mesh.nFy)[fyp]])
|
||||
sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm],
|
||||
np.arange(mesh.nFz)[fzp]])
|
||||
|
||||
xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
|
||||
xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
|
||||
@@ -165,6 +177,7 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
|
||||
|
||||
return xBC, yBC
|
||||
|
||||
|
||||
class Test1D_InhomogeneousMixed(Tests.OrderTest):
|
||||
name = "1D - Mixed"
|
||||
meshTypes = MESHTYPES
|
||||
@@ -173,11 +186,14 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest):
|
||||
meshSizes = [4, 8, 16, 32]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
phi_fun = lambda x: np.cos(np.pi*x)
|
||||
j_fun = lambda x: np.pi*np.sin(np.pi*x)
|
||||
phi_deriv = lambda x: -j_fun(x)
|
||||
q_fun = lambda x: (np.pi**2)*np.cos(np.pi*x)
|
||||
# Test function
|
||||
def phi_fun(x): return np.cos(np.pi*x)
|
||||
|
||||
def j_fun(x): return np.pi*np.sin(np.pi*x)
|
||||
|
||||
def phi_deriv(x): return -j_fun(x)
|
||||
|
||||
def q_fun(x): return (np.pi**2)*np.cos(np.pi*x)
|
||||
|
||||
xc_ana = phi_fun(self.M.gridCC)
|
||||
q_ana = q_fun(self.M.gridCC)
|
||||
@@ -194,12 +210,11 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest):
|
||||
beta = np.r_[beta_xm, beta_xp]
|
||||
vecN = self.M.vectorNx
|
||||
vecC = self.M.vectorCCx
|
||||
phi_bc = phi_fun(vecN[[0,-1]])
|
||||
phi_deriv_bc = phi_deriv(vecN[[0,-1]])
|
||||
phi_bc = phi_fun(vecN[[0, -1]])
|
||||
phi_deriv_bc = phi_deriv(vecN[[0, -1]])
|
||||
gamma = alpha*phi_bc + beta*phi_deriv_bc
|
||||
x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
|
||||
|
||||
|
||||
sigma = np.ones(self.M.nC)
|
||||
Mfrho = self.M.getFaceInnerProduct(1./sigma)
|
||||
MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
|
||||
@@ -214,7 +229,7 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest):
|
||||
A = Div*MfrhoI*G
|
||||
|
||||
if self.myTest == 'xc':
|
||||
#TODO: fix the null space
|
||||
# TODO: fix the null space
|
||||
Ainv = Solver(A)
|
||||
xc = Ainv*rhs
|
||||
err = np.linalg.norm((xc-xc_ana), np.inf)
|
||||
@@ -222,13 +237,13 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest):
|
||||
NotImplementedError
|
||||
return err
|
||||
|
||||
|
||||
def test_order(self):
|
||||
print "==== Testing Mixed boudary conduction for CC-problem ===="
|
||||
self.name = "1D"
|
||||
self.myTest = 'xc'
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class Test2D_InhomogeneousMixed(Tests.OrderTest):
|
||||
name = "2D - Mixed"
|
||||
meshTypes = MESHTYPES
|
||||
@@ -237,40 +252,59 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest):
|
||||
meshSizes = [4, 8, 16, 32]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
|
||||
j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
|
||||
j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
|
||||
phideriv_funX = lambda x: -j_funX(x)
|
||||
phideriv_funY = lambda x: -j_funY(x)
|
||||
q_fun = lambda x: +2*(np.pi**2)*phi_fun(x)
|
||||
# Test function
|
||||
def phi_fun(x):
|
||||
return np.cos(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1])
|
||||
|
||||
def j_funX(x):
|
||||
return +np.pi*np.sin(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1])
|
||||
|
||||
def j_funY(x):
|
||||
return +np.pi*np.cos(np.pi*x[:, 0])*np.sin(np.pi*x[:, 1])
|
||||
|
||||
def phideriv_funX(x):
|
||||
return -j_funX(x)
|
||||
|
||||
def phideriv_funY(x):
|
||||
return -j_funY(x)
|
||||
|
||||
def q_fun(x):
|
||||
return +2*(np.pi**2)*phi_fun(x)
|
||||
|
||||
xc_ana = phi_fun(self.M.gridCC)
|
||||
q_ana = q_fun(self.M.gridCC)
|
||||
jX_ana = j_funX(self.M.gridFx)
|
||||
jY_ana = j_funY(self.M.gridFy)
|
||||
j_ana = np.r_[jX_ana,jY_ana]
|
||||
j_ana = np.r_[jX_ana, jY_ana]
|
||||
|
||||
# Get boundary locations
|
||||
fxm,fxp,fym,fyp = self.M.faceBoundaryInd
|
||||
gBFxm = self.M.gridFx[fxm,:]
|
||||
gBFxp = self.M.gridFx[fxp,:]
|
||||
gBFym = self.M.gridFy[fym,:]
|
||||
gBFyp = self.M.gridFy[fyp,:]
|
||||
fxm, fxp, fym, fyp = self.M.faceBoundaryInd
|
||||
gBFxm = self.M.gridFx[fxm, :]
|
||||
gBFxp = self.M.gridFx[fxp, :]
|
||||
gBFym = self.M.gridFy[fym, :]
|
||||
gBFyp = self.M.gridFy[fyp, :]
|
||||
|
||||
# Setup Mixed B.C (alpha, beta, gamma)
|
||||
alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
|
||||
beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
|
||||
alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
|
||||
beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
|
||||
alpha_xm = np.ones_like(gBFxm[:, 0])
|
||||
alpha_xp = np.ones_like(gBFxp[:, 0])
|
||||
beta_xm = np.ones_like(gBFxm[:, 0])
|
||||
beta_xp = np.ones_like(gBFxp[:, 0])
|
||||
alpha_ym = np.ones_like(gBFym[:, 1])
|
||||
alpha_yp = np.ones_like(gBFyp[:, 1])
|
||||
beta_ym = np.ones_like(gBFym[:, 1])
|
||||
beta_yp = np.ones_like(gBFyp[:, 1])
|
||||
|
||||
phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
|
||||
phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
|
||||
|
||||
phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
|
||||
phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
|
||||
phiderivX_bc_xm = phideriv_funX(gBFxm)
|
||||
phiderivX_bc_xp = phideriv_funX(gBFxp)
|
||||
phiderivY_bc_ym = phideriv_funY(gBFym)
|
||||
phiderivY_bc_yp = phideriv_funY(gBFyp)
|
||||
|
||||
def gamma_fun(alpha, beta, phi, phi_deriv):
|
||||
return alpha*phi + beta*phi_deriv
|
||||
|
||||
gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
|
||||
gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
|
||||
gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
|
||||
gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
|
||||
@@ -282,7 +316,6 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest):
|
||||
|
||||
x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
|
||||
|
||||
|
||||
sigma = np.ones(self.M.nC)
|
||||
Mfrho = self.M.getFaceInnerProduct(1./sigma)
|
||||
MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
|
||||
@@ -303,13 +336,13 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest):
|
||||
NotImplementedError
|
||||
return err
|
||||
|
||||
|
||||
def test_order(self):
|
||||
print "==== Testing Mixed boudary conduction for CC-problem ===="
|
||||
self.name = "2D"
|
||||
self.myTest = 'xc'
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class Test3D_InhomogeneousMixed(Tests.OrderTest):
|
||||
name = "3D - Mixed"
|
||||
meshTypes = MESHTYPES
|
||||
@@ -318,51 +351,74 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
|
||||
meshSizes = [4, 8, 16]
|
||||
|
||||
def getError(self):
|
||||
#Test function
|
||||
phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
|
||||
j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
|
||||
j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
|
||||
j_funZ = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.sin(np.pi*x[:,2])
|
||||
# Test function
|
||||
def phi_fun(x):
|
||||
return (np.cos(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) *
|
||||
np.cos(np.pi*x[:, 2]))
|
||||
|
||||
phideriv_funX = lambda x: -j_funX(x)
|
||||
phideriv_funY = lambda x: -j_funY(x)
|
||||
phideriv_funZ = lambda x: -j_funZ(x)
|
||||
def j_funX(x):
|
||||
return (np.pi*np.sin(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) *
|
||||
np.cos(np.pi*x[:, 2]))
|
||||
|
||||
q_fun = lambda x: 3*(np.pi**2)*phi_fun(x)
|
||||
def j_funY(x):
|
||||
return (np.pi*np.cos(np.pi*x[:, 0])*np.sin(np.pi*x[:, 1]) *
|
||||
np.cos(np.pi*x[:, 2]))
|
||||
|
||||
def j_funZ(x):
|
||||
return (np.pi*np.cos(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) *
|
||||
np.sin(np.pi*x[:, 2]))
|
||||
|
||||
def phideriv_funX(x): return -j_funX(x)
|
||||
|
||||
def phideriv_funY(x): return -j_funY(x)
|
||||
|
||||
def phideriv_funZ(x): return -j_funZ(x)
|
||||
|
||||
def q_fun(x): return 3*(np.pi**2)*phi_fun(x)
|
||||
|
||||
xc_ana = phi_fun(self.M.gridCC)
|
||||
q_ana = q_fun(self.M.gridCC)
|
||||
jX_ana = j_funX(self.M.gridFx)
|
||||
jY_ana = j_funY(self.M.gridFy)
|
||||
j_ana = np.r_[jX_ana,jY_ana,jY_ana]
|
||||
j_ana = np.r_[jX_ana, jY_ana, jY_ana]
|
||||
|
||||
# Get boundary locations
|
||||
fxm,fxp,fym,fyp,fzm,fzp = self.M.faceBoundaryInd
|
||||
gBFxm = self.M.gridFx[fxm,:]
|
||||
gBFxp = self.M.gridFx[fxp,:]
|
||||
gBFym = self.M.gridFy[fym,:]
|
||||
gBFyp = self.M.gridFy[fyp,:]
|
||||
gBFzm = self.M.gridFz[fzm,:]
|
||||
gBFzp = self.M.gridFz[fzp,:]
|
||||
fxm, fxp, fym, fyp, fzm, fzp = self.M.faceBoundaryInd
|
||||
gBFxm = self.M.gridFx[fxm, :]
|
||||
gBFxp = self.M.gridFx[fxp, :]
|
||||
gBFym = self.M.gridFy[fym, :]
|
||||
gBFyp = self.M.gridFy[fyp, :]
|
||||
gBFzm = self.M.gridFz[fzm, :]
|
||||
gBFzp = self.M.gridFz[fzp, :]
|
||||
|
||||
# Setup Mixed B.C (alpha, beta, gamma)
|
||||
alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
|
||||
beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
|
||||
alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
|
||||
beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
|
||||
alpha_zm, alpha_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
|
||||
beta_zm, beta_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
|
||||
|
||||
alpha_xm = np.ones_like(gBFxm[:, 0])
|
||||
alpha_xp = np.ones_like(gBFxp[:, 0])
|
||||
beta_xm = np.ones_like(gBFxm[:, 0])
|
||||
beta_xp = np.ones_like(gBFxp[:, 0])
|
||||
alpha_ym = np.ones_like(gBFym[:, 1])
|
||||
alpha_yp = np.ones_like(gBFyp[:, 1])
|
||||
beta_ym = np.ones_like(gBFym[:, 1])
|
||||
beta_yp = np.ones_like(gBFyp[:, 1])
|
||||
alpha_zm = np.ones_like(gBFzm[:, 2])
|
||||
alpha_zp = np.ones_like(gBFzp[:, 2])
|
||||
beta_zm = np.ones_like(gBFzm[:, 2])
|
||||
beta_zp = np.ones_like(gBFzp[:, 2])
|
||||
|
||||
phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
|
||||
phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
|
||||
phi_bc_zm, phi_bc_zp = phi_fun(gBFzm), phi_fun(gBFzp)
|
||||
|
||||
phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
|
||||
phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
|
||||
phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funZ(gBFzm), phideriv_funZ(gBFzp)
|
||||
phiderivX_bc_xm = phideriv_funX(gBFxm)
|
||||
phiderivX_bc_xp = phideriv_funX(gBFxp)
|
||||
phiderivY_bc_ym = phideriv_funY(gBFym)
|
||||
phiderivY_bc_yp = phideriv_funY(gBFyp)
|
||||
phiderivY_bc_zm = phideriv_funZ(gBFzm)
|
||||
phiderivY_bc_zp = phideriv_funZ(gBFzp)
|
||||
|
||||
def gamma_fun(alpha, beta, phi, phi_deriv):
|
||||
return alpha*phi + beta*phi_deriv
|
||||
|
||||
gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
|
||||
gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
|
||||
gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
|
||||
gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
|
||||
@@ -376,7 +432,6 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
|
||||
|
||||
x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
|
||||
|
||||
|
||||
sigma = np.ones(self.M.nC)
|
||||
Mfrho = self.M.getFaceInnerProduct(1./sigma)
|
||||
MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
|
||||
@@ -390,7 +445,7 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
|
||||
A = Div*MfrhoI*G
|
||||
|
||||
if self.myTest == 'xc':
|
||||
#TODO: fix the null space
|
||||
# TODO: fix the null space
|
||||
Ainv = Solver(A)
|
||||
xc = Ainv*rhs
|
||||
err = np.linalg.norm((xc-xc_ana), np.inf)
|
||||
@@ -398,14 +453,11 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
|
||||
NotImplementedError
|
||||
return err
|
||||
|
||||
|
||||
def test_order(self):
|
||||
print "==== Testing Mixed boudary conduction for CC-problem ===="
|
||||
self.name = "3D"
|
||||
self.myTest = 'xc'
|
||||
self.orderTest()
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
@@ -2,8 +2,6 @@ import numpy as np
|
||||
import unittest
|
||||
from SimPEG.Utils import mkvc
|
||||
from SimPEG import Mesh, Tests
|
||||
import unittest
|
||||
|
||||
|
||||
MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh']
|
||||
TOLERANCES = [0.9, 0.5, 0.5]
|
||||
|
||||
@@ -52,7 +52,7 @@ class Tests(unittest.TestCase):
|
||||
assert o <= -1
|
||||
assert not (o > -1)
|
||||
assert o >= -1
|
||||
assert -(-o)*o == -o
|
||||
assert -1.*(-o)*o == -o
|
||||
o = Identity()
|
||||
assert +o == o
|
||||
assert -o == -o
|
||||
|
||||
Reference in New Issue
Block a user