From a950ba8dd1068f7107be342c11c210f8e8ef88da Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Mon, 18 Nov 2013 10:54:54 -0800 Subject: [PATCH 1/2] Bug fix in emSources. Beginning to change Cyl1DMesh to be more comparable with TensorMesh/BaseMesh. --- SimPEG/mesh/Cyl1DMesh.py | 22 +++++++++---------- .../utils/Geophysics/emSources/emSources.py | 1 + 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/SimPEG/mesh/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py index fa2bb729..28e2e1a9 100644 --- a/SimPEG/mesh/Cyl1DMesh.py +++ b/SimPEG/mesh/Cyl1DMesh.py @@ -62,11 +62,11 @@ class Cyl1DMesh(object): # Counting #################################################### - def nCr(): + def nCx(): doc = "Number of cells in the radial direction" fget = lambda self: self.hr.size return locals() - nCr = property(**nCr()) + nCx = property(**nCx()) def nCz(): doc = "Number of cells in the z direction" @@ -76,7 +76,7 @@ class Cyl1DMesh(object): def nC(): doc = "Total number of cells" - fget = lambda self: self.nCr * self.nCz + fget = lambda self: self.nCx * self.nCz return locals() nC = property(**nC()) @@ -106,7 +106,7 @@ class Cyl1DMesh(object): def nFz(): doc = "Number of z faces" - fget = lambda self: self.nNz * self.nCr + fget = lambda self: self.nNz * self.nCx return locals() nFz = property(**nFz()) @@ -242,12 +242,12 @@ class Cyl1DMesh(object): def fget(self): if self._edgeCurl is None: #1D Difference matricies - dr = sp.spdiags((np.ones((self.nCr+1, 1))*[-1, 1]).T, [-1,0], self.nCr, self.nCr, format="csr") + dr = sp.spdiags((np.ones((self.nCx+1, 1))*[-1, 1]).T, [-1,0], self.nCx, self.nCx, format="csr") dz = sp.spdiags((np.ones((self.nCz+1, 1))*[-1, 1]).T, [0,1], self.nCz, self.nCz+1, format="csr") #2D Difference matricies Dr = sp.kron(sp.eye(self.nNz), dr) - Dz = -sp.kron(dz, sp.eye(self.nCr)) #Not sure about this negative + Dz = -sp.kron(dz, sp.eye(self.nCx)) #Not sure about this negative #Edge curl operator self._edgeCurl = sp.diags(1/self.area,0)*sp.vstack((Dz, Dr))*sp.diags(self.edge,0) @@ -261,7 +261,7 @@ class Cyl1DMesh(object): def fget(self): if self._aveE2CC is None: az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr') - ar = sp.spdiags(0.5*np.ones((2, self.nCr)), [0, 1], self.nCr, self.nCr, format='csr') + ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr') ar[0,0] = 1 self._aveE2CC = sp.kron(az, ar).T return self._aveE2CC @@ -274,10 +274,10 @@ class Cyl1DMesh(object): def fget(self): if self._aveF2CC is None: az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr') - ar = sp.spdiags(0.5*np.ones((2, self.nCr)), [0, 1], self.nCr, self.nCr, format='csr') + ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr') ar[0,0] = 1 Afr = sp.kron(sp.eye(self.nCz),ar) - Afz = sp.kron(az,sp.eye(self.nCr)) + Afz = sp.kron(az,sp.eye(self.nCx)) self._aveF2CC = sp.vstack((Afr,Afz)).T return self._aveF2CC return locals() @@ -311,7 +311,7 @@ class Cyl1DMesh(object): elif type(materialProp) is float: materialProp = np.ones(self.nC)*materialProp elif materialProp.shape == (self.nCz,): - materialProp = materialProp.repeat(self.nCr) + materialProp = materialProp.repeat(self.nCx) materialProp = mkvc(materialProp) assert materialProp.shape == (self.nC,), "materialProp incorrect shape" @@ -383,7 +383,7 @@ class Cyl1DMesh(object): dFz = np.sum(dFz**2, axis=1) indBL = np.argmin(dFz) # Face below and to the left - indAL = indBL + self.nCr # Face above and to the left + indAL = indBL + self.nCx # Face above and to the left zF_BL = self.gridFz[indBL,:] zF_AL = self.gridFz[indAL,:] diff --git a/SimPEG/utils/Geophysics/emSources/emSources.py b/SimPEG/utils/Geophysics/emSources/emSources.py index 16c5a1fb..db492a94 100644 --- a/SimPEG/utils/Geophysics/emSources/emSources.py +++ b/SimPEG/utils/Geophysics/emSources/emSources.py @@ -25,6 +25,7 @@ def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0. txLoc = np.atleast_2d(txLoc) obsLoc = np.atleast_2d(obsLoc) + dipoleMoment = np.atleast_2d(dipoleMoment) nEdges = obsLoc.shape[0] nTx = txLoc.shape[0] From 8992f195a4008d96bf8d5a5a57559c4e7f6a4224 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Mon, 18 Nov 2013 11:06:35 -0800 Subject: [PATCH 2/2] Added nCv property to BaseMesh and Cyl1DMesh classes. --- SimPEG/mesh/BaseMesh.py | 11 +++++++++++ SimPEG/mesh/Cyl1DMesh.py | 6 ++++++ 2 files changed, 17 insertions(+) diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py index 2647cf44..ff804a83 100644 --- a/SimPEG/mesh/BaseMesh.py +++ b/SimPEG/mesh/BaseMesh.py @@ -228,6 +228,17 @@ class BaseMesh(object): return locals() nC = property(**nC()) + def nCv(): + doc = """ + Total number of cells in each direction + + :rtype: numpy.array (dim, ) + :return: [nCx, nCy, nCz] + """ + fget = lambda self: np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None]) + return locals() + nCv = property(**nCv()) + def nNx(): doc = """ Number of nodes in the x-direction diff --git a/SimPEG/mesh/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py index 28e2e1a9..e22e12b9 100644 --- a/SimPEG/mesh/Cyl1DMesh.py +++ b/SimPEG/mesh/Cyl1DMesh.py @@ -80,6 +80,12 @@ class Cyl1DMesh(object): return locals() nC = property(**nC()) + def nCv(): + doc = "Total number of cells in each direction" + fget = lambda self: np.array([self.nCx, self.nCz]) + return locals() + nCv = property(**nCv()) + def nNr(): doc = "Number of nodes in the radial direction" fget = lambda self: self.hr.size