From f8259456d03f2db29279730b432ae83c25f2d16c Mon Sep 17 00:00:00 2001 From: lheagy Date: Thu, 5 Feb 2015 12:44:46 -0800 Subject: [PATCH 1/4] start of cyl mesh anisotropy --- SimPEG/Maps.py | 24 ++++++++++++++++++++++++ SimPEG/Mesh/CylMesh.py | 25 +++++++++++++++++++++++++ 2 files changed, 49 insertions(+) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 4ab1d855..e37d435a 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -640,3 +640,27 @@ class CircleMap(IdentityMap): g4 = a*(-Y + y)*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1)*np.sqrt((X - x)**2 + (Y - y)**2)) g5 = -a*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1)) return np.c_[g1,g2,g3,g4,g5] + +class AnisotropyMap(IdentityMap): + + def __init__(self, mesh, **kwargs): + Utils.setKwargs(self, **kwargs) + self.mesh = mesh + + @property + def nP(self): + """ + :rtype: int + :return: number of parameters in the model + """ + return self.mesh.nC*2 + + @property + def shape(self): + """ + The default shape is (mesh.nC, nP). + + :rtype: (int,int) + :return: shape of the operator as a tuple + """ + return (self.mesh.nC*2, self.nP) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 93cd5451..d3a40dec 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -278,6 +278,19 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView): # kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") return self._aveE2CC + @property + def aveE2CCV(self): + "Construct the averaging operator on cell edges to cell centers." + if getattr(self, '_aveE2CCV', None) is None: + # The number of cell centers in each direction + n = self.vnC + if self.isSymmetric: + self._aveE2CCV = sp.block_diag((sp.kron(av(n[1]), speye(n[0])), + sp.kron(speye(n[1]), av(n[0]))), format="csr") + else: + raise NotImplementedError('wrapping in the averaging is not yet implemented') + return self._aveE2CCV + @property def aveF2CC(self): @@ -295,6 +308,18 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView): # kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") return self._aveF2CC + @property + def aveF2CCV(self): + "Construct the averaging operator on cell faces to cell centers." + if getattr(self, '_aveF2CCV', None) is None: + n = self.vnC + if self.isSymmetric: + self._aveF2CCV = sp.block_diag((sp.kron(speye(n[1]), av(n[0])), + sp.kron(av(n[1]), speye(n[0]))), format="csr") + else: + raise NotImplementedError('wrapping in the averaging is not yet implemented') + return self._aveF2CCV + if __name__ == '__main__': From 7069ca745134f484f91e0f8e8be24d6f5288700f Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 5 Feb 2015 13:51:30 -0800 Subject: [PATCH 2/4] updates to cyl averaging. --- SimPEG/Mesh/CylMesh.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index d3a40dec..8d2b1ae8 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -11,6 +11,10 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView): """ CylMesh is a mesh class for cylindrical problems + .. note:: + + for a cylindrically symmetric mesh use [hx, 1, hz] + :: cs, nc, npad = 20., 30, 8 @@ -25,8 +29,10 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView): def __init__(self, h, x0=None): BaseTensorMesh.__init__(self, h, x0) - assert self.dim == 3, "dim of mesh must equal 3, for a cylindrically symmetric mesh use [hx, 1, hz]" assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi" + if self.dim == 2: + print 'Warning, a disk mesh has not been tested thoroughly.' + @property def isSymmetric(self): @@ -285,8 +291,7 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView): # The number of cell centers in each direction n = self.vnC if self.isSymmetric: - self._aveE2CCV = sp.block_diag((sp.kron(av(n[1]), speye(n[0])), - sp.kron(speye(n[1]), av(n[0]))), format="csr") + return self.aveE2CC else: raise NotImplementedError('wrapping in the averaging is not yet implemented') return self._aveE2CCV @@ -314,8 +319,10 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView): if getattr(self, '_aveF2CCV', None) is None: n = self.vnC if self.isSymmetric: - self._aveF2CCV = sp.block_diag((sp.kron(speye(n[1]), av(n[0])), - sp.kron(av(n[1]), speye(n[0]))), format="csr") + avR = av(n[0])[:,1:] + avR[0,0] = 1. + self._aveF2CCV = sp.block_diag((sp.kron(speye(n[2]), avR), + sp.kron(av(n[2]), speye(n[0]))), format="csr") else: raise NotImplementedError('wrapping in the averaging is not yet implemented') return self._aveF2CCV From 708de19eb735ed247db06d2990364a3976b93919 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 5 Feb 2015 14:04:23 -0800 Subject: [PATCH 3/4] update the miniconda build. --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index dac2c6af..e31b94e3 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,7 +4,7 @@ python: # Setup anaconda before_install: - - if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.3.0-Linux-x86_64.sh -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.3.0-Linux-x86_64.sh -O miniconda.sh; fi + - if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh -O miniconda.sh; fi - chmod +x miniconda.sh - ./miniconda.sh -b - export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH From 28e0db4963aa059200998b4250a4245512260370 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 5 Feb 2015 14:22:20 -0800 Subject: [PATCH 4/4] remove anisotropy map --- SimPEG/Maps.py | 30 +++--------------------------- 1 file changed, 3 insertions(+), 27 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index e37d435a..0b47de66 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -229,16 +229,16 @@ class LogMap(IdentityMap): ..math:: - p = \\log(m) + p = \\log(m) - and + and ..math:: m = \\exp(p) NOTE: If you have a model which is log conductivity (ie. \\(m = \\log(\\sigma)\\)), - you should be using an ExpMap + you should be using an ExpMap """ @@ -640,27 +640,3 @@ class CircleMap(IdentityMap): g4 = a*(-Y + y)*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1)*np.sqrt((X - x)**2 + (Y - y)**2)) g5 = -a*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1)) return np.c_[g1,g2,g3,g4,g5] - -class AnisotropyMap(IdentityMap): - - def __init__(self, mesh, **kwargs): - Utils.setKwargs(self, **kwargs) - self.mesh = mesh - - @property - def nP(self): - """ - :rtype: int - :return: number of parameters in the model - """ - return self.mesh.nC*2 - - @property - def shape(self): - """ - The default shape is (mesh.nC, nP). - - :rtype: (int,int) - :return: shape of the operator as a tuple - """ - return (self.mesh.nC*2, self.nP)