diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py index aa487415..64a4f8b3 100644 --- a/SimPEG/Mesh/BaseMesh.py +++ b/SimPEG/Mesh/BaseMesh.py @@ -98,8 +98,7 @@ class BaseMesh(object): :rtype: int :return: nEy """ - if self.dim < 2: return None - return (self._n + np.r_[1,0,1][:self.dim]).prod() + return None if self.dim < 2 else (self._n + np.r_[1,0,1][:self.dim]).prod() @property def nEz(self): @@ -109,8 +108,7 @@ class BaseMesh(object): :rtype: int :return: nEz """ - if self.dim < 3: return None - return (self._n + np.r_[1,1,0][:self.dim]).prod() + return None if self.dim < 3 else (self._n + np.r_[1,1,0][:self.dim]).prod() @property def vnE(self): @@ -157,8 +155,7 @@ class BaseMesh(object): :rtype: int :return: nFy """ - if self.dim < 2: return None - return (self._n + np.r_[0,1,0][:self.dim]).prod() + return None if self.dim < 2 else (self._n + np.r_[0,1,0][:self.dim]).prod() @property def nFz(self): @@ -168,8 +165,7 @@ class BaseMesh(object): :rtype: int :return: nFz """ - if self.dim < 3: return None - return (self._n + np.r_[0,0,1][:self.dim]).prod() + return None if self.dim < 3 else (self._n + np.r_[0,0,1][:self.dim]).prod() @property def vnF(self): @@ -316,7 +312,7 @@ class BaseRectangularMesh(BaseMesh): @property def nNy(self): """ - Number of noes in the y-direction + Number of nodes in the y-direction :rtype: int :return: nNy or None if dim < 2 diff --git a/SimPEG/Mesh/Cyl1DMesh.py b/SimPEG/Mesh/Cyl1DMesh.py deleted file mode 100644 index f27d3d0f..00000000 --- a/SimPEG/Mesh/Cyl1DMesh.py +++ /dev/null @@ -1,476 +0,0 @@ -import numpy as np -import scipy.sparse as sp -from scipy.constants import pi -from SimPEG.Utils import mkvc, ndgrid, sdiag - -class Cyl1DMesh(object): - """ - Cyl1DMesh is a mesh class for cylindrically symmetric 1D problems - """ - - _meshType = 'CYL1D' - - def __init__(self, h, z0=None): - assert len(h) == 2, "len(h) must equal 2" - if z0 is not None: - assert z0.size == 1, "z0.size must equal 1" - - for i, h_i in enumerate(h): - assert type(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) - - # Ensure h contains 1D vectors - self._h = [mkvc(x.astype(float)) for x in h] - - if z0 is None: - z0 = 0 - self._z0 = z0 - - #################################################### - # Mesh properties - #################################################### - - def h(): - doc = "list containing the width of each cell" - def fget(self): - return self._h - return locals() - h = property(**h()) - - @property - def dim(self): return 2 - - def z0(): - doc = "The z-origin" - def fget(self): - return self._z0 - return locals() - z0 = property(**z0()) - - def hr(): - doc = "Width of the cells in the r direction" - def fget(self): - return self._h[0] - return locals() - hr = property(**hr()) - - def hz(): - doc = "Width of the cells in the z direction" - def fget(self): - return self._h[1] - return locals() - hz = property(**hz()) - - #################################################### - # Counting - #################################################### - - def nCx(): - doc = "Number of cells in the radial direction" - fget = lambda self: self.hr.size - return locals() - nCx = property(**nCx()) - - def nCz(): - doc = "Number of cells in the z direction" - fget = lambda self: self.hz.size - return locals() - nCz = property(**nCz()) - - def nC(): - doc = "Total number of cells" - fget = lambda self: self.nCx * self.nCz - return locals() - nC = property(**nC()) - - def vnC(): - doc = "Total number of cells in each direction" - fget = lambda self: np.array([self.nCx, self.nCz]) - return locals() - vnC = property(**vnC()) - - def nNr(): - doc = "Number of nodes in the radial direction" - fget = lambda self: self.hr.size - return locals() - nNr = property(**nNr()) - - def nNz(): - doc = "Number of nodes in the radial direction" - fget = lambda self: self.hz.size + 1 - return locals() - nNz = property(**nNz()) - - def nN(): - doc = "Total number of nodes" - fget = lambda self: self.nNr * self.nNz - return locals() - nN = property(**nN()) - - def nFr(): - doc = "Number of r faces" - fget = lambda self: self.nNr * self.nCz - return locals() - nFr = property(**nFr()) - - def vnFz(): - doc = "Number of z faces" - fget = lambda self: self.nNz * self.nCx - return locals() - vnFz = property(**vnFz()) - - def vnF(): - doc = "Total number of faces in each direction" - fget = lambda self: np.array([self.nFr, self.vnFz]) - return locals() - vnF = property(**vnF()) - - def nF(): - doc = "Total number of faces" - fget = lambda self: self.nFr + self.vnFz - return locals() - nF = property(**nF()) - - def nE(): - doc = "Number of edges" - fget = lambda self: self.nN - return locals() - nE = property(**nE()) - - #################################################### - # Vectors & Grids - #################################################### - - def vectorNr(): - doc = "Nodal grid vector (1D) in the r direction" - fget = lambda self: self.hr.cumsum() - return locals() - vectorNr = property(**vectorNr()) - - def vectorNz(): - doc = "Nodal grid vector (1D) in the z direction" - fget = lambda self: np.r_[0, self.hz.cumsum()] + self._z0 - return locals() - vectorNz = property(**vectorNz()) - - def vectorCCr(): - doc = "Cell centered grid vector (1D) in the r direction" - fget = lambda self: np.r_[0, self.hr.cumsum()[1:] - self.hr[1:]/2] - return locals() - vectorCCr = property(**vectorCCr()) - - def vectorCCz(): - doc = "Cell centered grid vector (1D) in the z direction" - fget = lambda self: self.hz.cumsum() - self.hz/2 + self._z0 - return locals() - vectorCCz = property(**vectorCCz()) - - def gridCC(): - doc = "Cell-centered grid" - def fget(self): - if self._gridCC is None: - self._gridCC = ndgrid([self.vectorCCr, self.vectorCCz]) - return self._gridCC - return locals() - _gridCC = None - gridCC = property(**gridCC()) - - def gridN(): - doc = "Nodal grid" - def fget(self): - if self._gridN is None: - self._gridN = ndgrid([self.vectorNr, self.vectorNz]) - return self._gridN - return locals() - _gridN = None - gridN = property(**gridN()) - - def gridFr(): - doc = "r face grid" - def fget(self): - if self._gridFr is None: - self._gridFr = ndgrid([self.vectorNr, self.vectorCCz]) - return self._gridFr - return locals() - _gridFr = None - gridFr = property(**gridFr()) - - def gridFz(): - doc = "z face grid" - def fget(self): - if self._gridFz is None: - self._gridFz = ndgrid([self.vectorCCr, self.vectorNz]) - return self._gridFz - return locals() - _gridFz = None - gridFz = property(**gridFz()) - - #################################################### - # Geometries - #################################################### - - def edge(): - doc = "Edge lengths" - def fget(self): - if self._edge is None: - self._edge = 2*pi*self.gridN[:,0] - return self._edge - return locals() - _edge = None - edge = property(**edge()) - - def area(): - doc = "Face areas" - def fget(self): - if self._area is None: - areaR = np.kron(self.hz, 2*pi*self.vectorNr) - areaZ = np.kron(np.ones_like(self.vectorNz),pi*(self.vectorNr**2 - np.r_[0, self.vectorNr[:-1]]**2)) - self._area = np.r_[areaR, areaZ] - return self._area - return locals() - _area = None - area = property(**area()) - - def vol(): - doc = "Volume of each cell" - def fget(self): - if self._vol is None: - az = pi*(self.vectorNr**2 - np.r_[0, self.vectorNr[:-1]]**2) - self._vol = np.kron(self.hz,az) - return self._vol - return locals() - _vol = None - vol = property(**vol()) - - #################################################### - # Operators - #################################################### - - def edgeCurl(): - doc = "The edgeCurl property." - def fget(self): - if self._edgeCurl is None: - #1D Difference matricies - 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.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) - return self._edgeCurl - return locals() - _edgeCurl = None - edgeCurl = property(**edgeCurl()) - - def aveE2CC(): - doc = "Averaging operator from cell edges to cell centres" - 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.nCx)), [0, 1], self.nCx, self.nCx, format='csr') - ar[0,0] = 1 - self._aveE2CC = sp.kron(az, ar).T - return self._aveE2CC - return locals() - _aveE2CC = None - aveE2CC = property(**aveE2CC()) - - def aveF2CC(): - doc = "Averaging operator from cell faces to cell centres" - 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.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.nCx)) - self._aveF2CC = sp.vstack((Afr,Afz)).T - return self._aveF2CC - return locals() - _aveF2CC = None - aveF2CC = property(**aveF2CC()) - - def getFaceMassDeriv(self): - Av = self.aveF2CC - return Av.T * sdiag(self.vol) - - def getEdgeMassDeriv(self): - Av = self.aveE2CC - return Av.T * sdiag(self.vol) - - - #################################################### - # Methods - #################################################### - - - def getMass(self, materialProp=None, loc='e'): - """ Produces mass matricies. - - :param None,float,numpy.ndarray materialProp: property to be averaged (see below) - :param str loc: Average to location: 'e'-edges, 'f'-faces - :rtype: scipy.sparse.csr.csr_matrix - :return: M, the mass matrix - - materialProp can be:: - - None -> takes materialProp = 1 (default) - float -> a constant value for entire domain - numpy.ndarray -> if materialProp.size == self.nC - 3D property model - if materialProp.size = self.nCz - 1D (layered eath) property model - """ - if materialProp is None: - materialProp = np.ones(self.nC) - elif type(materialProp) is float: - materialProp = np.ones(self.nC)*materialProp - elif materialProp.shape == (self.nCz,): - materialProp = materialProp.repeat(self.nCx) - materialProp = mkvc(materialProp) - assert materialProp.shape == (self.nC,), "materialProp incorrect shape" - - if loc=='e': - Av = self.aveE2CC - elif loc=='f': - Av = self.aveF2CC - else: - raise ValueError('Invalid loc') - - diag = Av.T * (self.vol * mkvc(materialProp)) - - return sdiag(diag) - - def getEdgeMass(self, materialProp=None): - """mass matrix for products of edge functions w'*M(materialProp)*e""" - return self.getMass(loc='e', materialProp=materialProp) - - def getFaceMass(self, materialProp=None): - """mass matrix for products of face functions w'*M(materialProp)*f""" - return self.getMass(loc='f', materialProp=materialProp) - - def getInterpolationMat(self, loc, locType='fz'): - """ Produces intrpolation matrix - - :param numpy.ndarray loc: Location of points to interpolate to - :param str locType: What to interpolate (see below) - :rtype: scipy.sparse.csr.csr_matrix - :return: M, the intrpolation matrix - - locType can be:: - - 'fz' -> z-component of field defined on faces - 'fr' -> r-component of field defined on faces - 'et' -> theta-component of field defined on edges - """ - - loc = np.atleast_2d(loc) - - assert np.all(loc[:,0]<=self.vectorNr.max()) & \ - np.all(loc[:,1]>=self.vectorNz.min()) & \ - np.all(loc[:,1]<=self.vectorNz.max()), \ - "Points outside of mesh" - - - if locType=='fz': - Q = sp.lil_matrix((loc.shape[0], self.nF), dtype=float) - - for i, iloc in enumerate(loc): - # Point is on a z-interface - if np.any(np.abs(self.vectorNz-iloc[1])<0.001): - dFz = self.gridFz-iloc #Distance to z faces - dFz[dFz[:,0]>0,:] = np.inf #Looking for next face to the left... - indL = np.argmin(np.sum(dFz**2, axis=1)) #Closest one - if self.gridFz[indL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation) - zFL = self.gridFz[indL,:] - zFLL = self.gridFz[indL-1,:] - Q[i, indL+self.nFr] = (iloc[0] - zFLL[0])/(zFL[0] - zFLL[0]) - Q[i, indL+self.nFr-1] = -(iloc[0] - zFL[0])/(zFL[0] - zFLL[0]) - else: - zFL = self.gridFz[indL,:] - zFR = self.gridFz[indL+1,:] - Q[i,indL+self.nFr] = (zFR[0] - iloc[0])/(zFR[0] - zFL[0]) - Q[i,indL+self.nFr+1] = (iloc[0] - zFL[0])/(zFR[0] - zFL[0]) - # Point is in a cell - else: - dFz = self.gridFz-iloc - dFz[dFz>0] = np.inf - dFz = np.sum(dFz**2, axis=1) - - indBL = np.argmin(dFz) # Face below and to the left - indAL = indBL + self.nCx # Face above and to the left - - zF_BL = self.gridFz[indBL,:] - zF_AL = self.gridFz[indAL,:] - - dzB = iloc[1] - zF_BL[1] # z-distance to face below - dzA = zF_AL[1] - iloc[1] # z-distance to face above - - if self.gridFz[indBL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation) - zF_BLL = self.gridFz[indBL-1,:] - zF_ALL = self.gridFz[indAL-1,:] - - DZ = zF_AL[1] - zF_BL[1] - DR = zF_AL[0] - zF_ALL[0] - - drL = iloc[0] - zF_AL[0] - drLL = iloc[0] - zF_ALL[0] - - Q[i, indBL+self.nFr-1] = -(1 - dzB/DZ)*(drL/DR) - Q[i, indBL+self.nFr] = (1 - dzB/DZ)*(drLL/DR) - Q[i, indAL+self.nFr-1] = -(dzB/DZ)*(drL/DR) - Q[i, indAL+self.nFr] = (dzB/DZ)*(drLL/DR) - else: - indBR = indBL+1 # Face below and to the right - indAR = indAL + 1 # Face above and to the right - zF_BR = self.gridFz[indBR,:] - - drL = iloc[0] - zF_BL[0] # r-distance to face on left - drR = zF_BR[0] - iloc[0] # r-distance to face on right - - drz = (drL + drR)*(dzB + dzA) - Q[i,indBL+self.nFr] = drR*dzA/drz - Q[i,indBR+self.nFr] = drL*dzA/drz - Q[i,indAL+self.nFr] = drR*dzB/drz - Q[i,indAR+self.nFr] = drL*dzB/drz - - elif locType=='fr': - raise NotImplementedError('locType==fr') - elif locType=='et': - raise NotImplementedError('locType==et') - else: - raise ValueError('Invalid locType') - return Q.tocsr() - - def getNearest(self, loc, locType): - """ Returns the index of the closest face or edge to a given location - - :param numpy.ndarray loc: Test point - :param str locType: Type of location desired (see below) - :rtype: int - :return: ind: - - locType can be:: - - 'fz' -> location of nearest z-face - 'fr' -> location of nearest r-face - 'et' -> location of nearest edge - """ - - if locType=='et': - dr = self.gridN[:,0] - loc[0] - dz = self.gridN[:,1] - loc[1] - elif locType=='fz': - dr = self.gridFz[:,0] - loc[0] - dz = self.gridFz[:,1] - loc[1] - elif locType=='fr': - dr = self.gridFr[:,0] - loc[0] - dz = self.gridFr[:,1] - loc[1] - else: - raise ValueError('Invalid locType') - R = np.sqrt(dr**2 + dz**2) - ind = np.argmin(R) - return ind diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py new file mode 100644 index 00000000..594c739c --- /dev/null +++ b/SimPEG/Mesh/CylMesh.py @@ -0,0 +1,337 @@ +import numpy as np +import scipy.sparse as sp +from scipy.constants import pi +from SimPEG.Utils import mkvc, ndgrid, sdiag +from TensorMesh import TensorMesh + +class CylMesh(TensorMesh): + """ + CylMesh is a mesh class for cylindrically problems + """ + + _meshType = 'CYL' + + def __init__(self, h, x0=None): + assert len(h) == 3, "len(h) must equal 3, for a cylindrically symmetric mesh use [hx, 1, hz]" + + if x0 is not None: + assert x0.size == 3, "x0.size must equal 1" + else: + x0 = np.r_[0, 0, 0] + + for i, h_i in enumerate(h): + if type(h_i) in [int, long, float]: + # This gives you something over the unit cylinder. + h_i = (2*np.pi if i==1 else 1.)*np.ones(int(h_i))/int(h_i) + assert type(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) + h[i] = h_i[:] # make a copy. + + assert h[1].sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi" + + TensorMesh.__init__(self, h, x0) + + + @property + def nNy(self): + """ + Number of nodes in the y-direction + + :rtype: int + :return: nNy + """ + return self.nCy + + @property + def nN(self): + """ + Total number of nodes + + :rtype: int + :return: nN + + .. plot:: + :include-source: + + from SimPEG import Mesh, np + Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(nodes=True,showIt=True) + """ + return (np.r_[self.nNx, self.nNy, self.nNz]).prod() + + @property + def nFy(self): + """ + Number of y-faces + + :rtype: int + :return: nFy + """ + return (self.vnC + np.r_[0,-1,0][:self.dim]).prod() + + @property + def vnF(self): + """Total number of faces in each direction""" + return np.array([self.nFr, self.vnFz]) + + @property + def nF(self): + """Total number of faces""" + return self.nFr + self.vnFz + + @property + def nE(self): + """Number of edges""" + return self.nN + + @property + def vectorNx(self): + """Nodal grid vector (1D) in the r direction""" + return self.hr.cumsum() + + @property + def edge(self): + """Edge lengths""" + if getattr(self, '_edge', None) is None: + self._edge = 2*pi*self.gridN[:,0] + return self._edge + + @property + def area(self): + """Face areas""" + if getattr(self, '_area', None) is None: + areaR = np.kron(self.hz, 2*pi*self.vectorNr) + areaZ = np.kron(np.ones_like(self.vectorNz),pi*(self.vectorNr**2 - np.r_[0, self.vectorNr[:-1]]**2)) + self._area = np.r_[areaR, areaZ] + return self._area + + @property + def vol(self): + """Volume of each cell""" + if getattr(self, '_vol', None) is None: + az = pi*(self.vectorNr**2 - np.r_[0, self.vectorNr[:-1]]**2) + self._vol = np.kron(self.hz,az) + return self._vol + + #################################################### + # Operators + #################################################### + + @property + def edgeCurl(self): + """The edgeCurl property.""" + if getattr(self, '_edgeCurl', None) is None: + #1D Difference matricies + 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.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) + return self._edgeCurl + + @property + def aveE2CC(self): + """Averaging operator from cell edges to cell centres""" + if getattr(self, '_aveE2CC', None) 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.nCx)), [0, 1], self.nCx, self.nCx, format='csr') + ar[0,0] = 1 + self._aveE2CC = sp.kron(az, ar).T + return self._aveE2CC + + @property + def aveF2CC(self): + """Averaging operator from cell faces to cell centres""" + if getattr(self, '_aveF2CC', None) 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.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.nCx)) + self._aveF2CC = sp.vstack((Afr,Afz)).T + return self._aveF2CC + + def getFaceMassDeriv(self): + Av = self.aveF2CC + return Av.T * sdiag(self.vol) + + def getEdgeMassDeriv(self): + Av = self.aveE2CC + return Av.T * sdiag(self.vol) + + + #################################################### + # Methods + #################################################### + + + def getMass(self, materialProp=None, loc='e'): + """ Produces mass matricies. + + :param None,float,numpy.ndarray materialProp: property to be averaged (see below) + :param str loc: Average to location: 'e'-edges, 'f'-faces + :rtype: scipy.sparse.csr.csr_matrix + :return: M, the mass matrix + + materialProp can be:: + + None -> takes materialProp = 1 (default) + float -> a constant value for entire domain + numpy.ndarray -> if materialProp.size == self.nC + 3D property model + if materialProp.size = self.nCz + 1D (layered eath) property model + """ + if materialProp is None: + materialProp = np.ones(self.nC) + elif type(materialProp) is float: + materialProp = np.ones(self.nC)*materialProp + elif materialProp.shape == (self.nCz,): + materialProp = materialProp.repeat(self.nCx) + materialProp = mkvc(materialProp) + assert materialProp.shape == (self.nC,), "materialProp incorrect shape" + + if loc=='e': + Av = self.aveE2CC + elif loc=='f': + Av = self.aveF2CC + else: + raise ValueError('Invalid loc') + + diag = Av.T * (self.vol * mkvc(materialProp)) + + return sdiag(diag) + + def getEdgeMass(self, materialProp=None): + """mass matrix for products of edge functions w'*M(materialProp)*e""" + return self.getMass(loc='e', materialProp=materialProp) + + def getFaceMass(self, materialProp=None): + """mass matrix for products of face functions w'*M(materialProp)*f""" + return self.getMass(loc='f', materialProp=materialProp) + + def getInterpolationMat(self, loc, locType='fz'): + """ Produces intrpolation matrix + + :param numpy.ndarray loc: Location of points to interpolate to + :param str locType: What to interpolate (see below) + :rtype: scipy.sparse.csr.csr_matrix + :return: M, the intrpolation matrix + + locType can be:: + + 'fz' -> z-component of field defined on faces + 'fr' -> r-component of field defined on faces + 'et' -> theta-component of field defined on edges + """ + + loc = np.atleast_2d(loc) + + assert np.all(loc[:,0]<=self.vectorNr.max()) & \ + np.all(loc[:,1]>=self.vectorNz.min()) & \ + np.all(loc[:,1]<=self.vectorNz.max()), \ + "Points outside of mesh" + + + if locType=='fz': + Q = sp.lil_matrix((loc.shape[0], self.nF), dtype=float) + + for i, iloc in enumerate(loc): + # Point is on a z-interface + if np.any(np.abs(self.vectorNz-iloc[1])<0.001): + dFz = self.gridFz-iloc #Distance to z faces + dFz[dFz[:,0]>0,:] = np.inf #Looking for next face to the left... + indL = np.argmin(np.sum(dFz**2, axis=1)) #Closest one + if self.gridFz[indL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation) + zFL = self.gridFz[indL,:] + zFLL = self.gridFz[indL-1,:] + Q[i, indL+self.nFr] = (iloc[0] - zFLL[0])/(zFL[0] - zFLL[0]) + Q[i, indL+self.nFr-1] = -(iloc[0] - zFL[0])/(zFL[0] - zFLL[0]) + else: + zFL = self.gridFz[indL,:] + zFR = self.gridFz[indL+1,:] + Q[i,indL+self.nFr] = (zFR[0] - iloc[0])/(zFR[0] - zFL[0]) + Q[i,indL+self.nFr+1] = (iloc[0] - zFL[0])/(zFR[0] - zFL[0]) + # Point is in a cell + else: + dFz = self.gridFz-iloc + dFz[dFz>0] = np.inf + dFz = np.sum(dFz**2, axis=1) + + indBL = np.argmin(dFz) # Face below and to the left + indAL = indBL + self.nCx # Face above and to the left + + zF_BL = self.gridFz[indBL,:] + zF_AL = self.gridFz[indAL,:] + + dzB = iloc[1] - zF_BL[1] # z-distance to face below + dzA = zF_AL[1] - iloc[1] # z-distance to face above + + if self.gridFz[indBL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation) + zF_BLL = self.gridFz[indBL-1,:] + zF_ALL = self.gridFz[indAL-1,:] + + DZ = zF_AL[1] - zF_BL[1] + DR = zF_AL[0] - zF_ALL[0] + + drL = iloc[0] - zF_AL[0] + drLL = iloc[0] - zF_ALL[0] + + Q[i, indBL+self.nFr-1] = -(1 - dzB/DZ)*(drL/DR) + Q[i, indBL+self.nFr] = (1 - dzB/DZ)*(drLL/DR) + Q[i, indAL+self.nFr-1] = -(dzB/DZ)*(drL/DR) + Q[i, indAL+self.nFr] = (dzB/DZ)*(drLL/DR) + else: + indBR = indBL+1 # Face below and to the right + indAR = indAL + 1 # Face above and to the right + zF_BR = self.gridFz[indBR,:] + + drL = iloc[0] - zF_BL[0] # r-distance to face on left + drR = zF_BR[0] - iloc[0] # r-distance to face on right + + drz = (drL + drR)*(dzB + dzA) + Q[i,indBL+self.nFr] = drR*dzA/drz + Q[i,indBR+self.nFr] = drL*dzA/drz + Q[i,indAL+self.nFr] = drR*dzB/drz + Q[i,indAR+self.nFr] = drL*dzB/drz + + elif locType=='fr': + raise NotImplementedError('locType==fr') + elif locType=='et': + raise NotImplementedError('locType==et') + else: + raise ValueError('Invalid locType') + return Q.tocsr() + + def getNearest(self, loc, locType): + """ Returns the index of the closest face or edge to a given location + + :param numpy.ndarray loc: Test point + :param str locType: Type of location desired (see below) + :rtype: int + :return: ind: + + locType can be:: + + 'fz' -> location of nearest z-face + 'fr' -> location of nearest r-face + 'et' -> location of nearest edge + """ + + if locType=='et': + dr = self.gridN[:,0] - loc[0] + dz = self.gridN[:,1] - loc[1] + elif locType=='fz': + dr = self.gridFz[:,0] - loc[0] + dz = self.gridFz[:,1] - loc[1] + elif locType=='fr': + dr = self.gridFr[:,0] - loc[0] + dz = self.gridFr[:,1] - loc[1] + else: + raise ValueError('Invalid locType') + R = np.sqrt(dr**2 + dz**2) + ind = np.argmin(R) + return ind diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 12253bb5..693359bf 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -105,226 +105,172 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): return outStr - def h(): - doc = "h is a list containing the cell widths of the tensor mesh in each dimension." - fget = lambda self: self._h - return locals() - h = property(**h()) + @property + def h(self): + """h is a list containing the cell widths of the tensor mesh in each dimension.""" + return self._h - def hx(): - doc = "Width of cells in the x direction" - fget = lambda self: self._h[0] - return locals() - hx = property(**hx()) + @property + def hx(self): + "Width of cells in the x direction" + return self._h[0] - def hy(): - doc = "Width of cells in the y direction" - fget = lambda self: None if self.dim < 2 else self._h[1] - return locals() - hy = property(**hy()) + @property + def hy(self): + "Width of cells in the y direction" + return None if self.dim < 2 else self._h[1] - def hz(): - doc = "Width of cells in the z direction" - fget = lambda self: None if self.dim < 3 else self._h[2] - return locals() - hz = property(**hz()) + @property + def hz(self): + "Width of cells in the z direction" + return None if self.dim < 3 else self._h[2] - def vectorNx(): - doc = "Nodal grid vector (1D) in the x direction." - fget = lambda self: np.r_[0., self.hx.cumsum()] + self.x0[0] - return locals() - vectorNx = property(**vectorNx()) + @property + def vectorNx(self): + """Nodal grid vector (1D) in the x direction.""" + return np.r_[0., self.hx.cumsum()] + self.x0[0] - def vectorNy(): - doc = "Nodal grid vector (1D) in the y direction." - fget = lambda self: None if self.dim < 2 else np.r_[0., self.hy.cumsum()] + self.x0[1] - return locals() - vectorNy = property(**vectorNy()) + @property + def vectorNy(self): + """Nodal grid vector (1D) in the y direction.""" + return None if self.dim < 2 else np.r_[0., self.hy.cumsum()] + self.x0[1] - def vectorNz(): - doc = "Nodal grid vector (1D) in the z direction." - fget = lambda self: None if self.dim < 3 else np.r_[0., self.hz.cumsum()] + self.x0[2] - return locals() - vectorNz = property(**vectorNz()) + @property + def vectorNz(self): + """Nodal grid vector (1D) in the z direction.""" + return None if self.dim < 3 else np.r_[0., self.hz.cumsum()] + self.x0[2] - def vectorCCx(): - doc = "Cell-centered grid vector (1D) in the x direction." - fget = lambda self: np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + self.x0[0] - return locals() - vectorCCx = property(**vectorCCx()) + @property + def vectorCCx(self): + """Cell-centered grid vector (1D) in the x direction.""" + return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + self.x0[0] - def vectorCCy(): - doc = "Cell-centered grid vector (1D) in the y direction." - fget = lambda self: None if self.dim < 2 else np.r_[0, self.hy[:-1].cumsum()] + self.hy*0.5 + self.x0[1] - return locals() - vectorCCy = property(**vectorCCy()) + @property + def vectorCCy(self): + """Cell-centered grid vector (1D) in the y direction.""" + return None if self.dim < 2 else np.r_[0, self.hy[:-1].cumsum()] + self.hy*0.5 + self.x0[1] - def vectorCCz(): - doc = "Cell-centered grid vector (1D) in the z direction." - fget = lambda self: None if self.dim < 3 else np.r_[0, self.hz[:-1].cumsum()] + self.hz*0.5 + self.x0[2] - return locals() - vectorCCz = property(**vectorCCz()) + @property + def vectorCCz(self): + """Cell-centered grid vector (1D) in the z direction.""" + return None if self.dim < 3 else np.r_[0, self.hz[:-1].cumsum()] + self.hz*0.5 + self.x0[2] - def gridCC(): - doc = "Cell-centered grid." + @property + def gridCC(self): + """Cell-centered grid.""" + if getattr(self, '_gridCC', None) is None: + self._gridCC = Utils.ndgrid(self.getTensor('CC')) + return self._gridCC - def fget(self): - if self._gridCC is None: - self._gridCC = Utils.ndgrid(self.getTensor('CC')) - return self._gridCC - return locals() - _gridCC = None # Store grid by default - gridCC = property(**gridCC()) + @property + def gridN(self): + """Nodal grid.""" + if getattr(self, '_gridN', None) is None: + self._gridN = Utils.ndgrid(self.getTensor('N')) + return self._gridN - def gridN(): - doc = "Nodal grid." + @property + def gridFx(self): + """Face staggered grid in the x direction.""" + if getattr(self, '_gridFx', None) is None: + self._gridFx = Utils.ndgrid(self.getTensor('Fx')) + return self._gridFx - def fget(self): - if self._gridN is None: - self._gridN = Utils.ndgrid(self.getTensor('N')) - return self._gridN - return locals() - _gridN = None # Store grid by default - gridN = property(**gridN()) + @property + def gridFy(self): + """Face staggered grid in the y direction.""" + if getattr(self, '_gridFy', None) is None and self.dim > 1: + self._gridFy = Utils.ndgrid(self.getTensor('Fy')) + return self._gridFy - def gridFx(): - doc = "Face staggered grid in the x direction." + @property + def gridFz(self): + """Face staggered grid in the z direction.""" + if getattr(self, '_gridFz', None) is None and self.dim > 2: + self._gridFz = Utils.ndgrid(self.getTensor('Fz')) + return self._gridFz - def fget(self): - if self._gridFx is None: - self._gridFx = Utils.ndgrid(self.getTensor('Fx')) - return self._gridFx - return locals() - _gridFx = None # Store grid by default - gridFx = property(**gridFx()) + @property + def gridEx(self): + """Edge staggered grid in the x direction.""" + if getattr(self, '_gridEx', None) is None: + self._gridEx = Utils.ndgrid(self.getTensor('Ex')) + return self._gridEx - def gridFy(): - doc = "Face staggered grid in the y direction." + @property + def gridEy(self): + """Edge staggered grid in the y direction.""" + if getattr(self, '_gridEy', None) is None and self.dim > 1: + self._gridEy = Utils.ndgrid(self.getTensor('Ey')) + return self._gridEy - def fget(self): - if self._gridFy is None and self.dim > 1: - self._gridFy = Utils.ndgrid(self.getTensor('Fy')) - return self._gridFy - return locals() - _gridFy = None # Store grid by default - gridFy = property(**gridFy()) - - def gridFz(): - doc = "Face staggered grid in the z direction." - - def fget(self): - if self._gridFz is None and self.dim > 2: - self._gridFz = Utils.ndgrid(self.getTensor('Fz')) - 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: - self._gridEx = Utils.ndgrid(self.getTensor('Ex')) - 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 and self.dim > 1: - self._gridEy = Utils.ndgrid(self.getTensor('Ey')) - 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 > 2: - self._gridEz = Utils.ndgrid(self.getTensor('Ez')) - return self._gridEz - return locals() - _gridEz = None # Store grid by default - gridEz = property(**gridEz()) + @property + def gridEz(self): + """Edge staggered grid in the z direction.""" + if getattr(self, '_gridEz', None) is None and self.dim > 2: + self._gridEz = Utils.ndgrid(self.getTensor('Ez')) + return self._gridEz # --------------- Geometries --------------------- - def vol(): - doc = "Construct cell volumes of the 3D model as 1d array." + @property + def vol(self): + """Construct cell volumes of the 3D model as 1d array.""" + if getattr(self, '_vol', None) is None: + vh = self.h + # Compute cell volumes + if self.dim == 1: + self._vol = Utils.mkvc(vh[0]) + elif self.dim == 2: + # Cell sizes in each direction + self._vol = Utils.mkvc(np.outer(vh[0], vh[1])) + elif self.dim == 3: + # Cell sizes in each direction + self._vol = Utils.mkvc(np.outer(Utils.mkvc(np.outer(vh[0], vh[1])), vh[2])) + return self._vol - def fget(self): - if(self._vol is None): - vh = self.h - # Compute cell volumes - if(self.dim == 1): - self._vol = Utils.mkvc(vh[0]) - elif(self.dim == 2): - # Cell sizes in each direction - self._vol = Utils.mkvc(np.outer(vh[0], vh[1])) - elif(self.dim == 3): - # Cell sizes in each direction - self._vol = Utils.mkvc(np.outer(Utils.mkvc(np.outer(vh[0], vh[1])), vh[2])) - return self._vol - return locals() - _vol = None - vol = property(**vol()) + @property + def area(self): + """Construct face areas of the 3D model as 1d array.""" + if getattr(self, '_area', None) is None: + # Ensure that we are working with column vectors + vh = self.h + # The number of cell centers in each direction + n = self.vnC + # Compute areas of cell faces + if(self.dim == 1): + self._area = np.ones(n[0]+1) + elif(self.dim == 2): + area1 = np.outer(np.ones(n[0]+1), vh[1]) + area2 = np.outer(vh[0], np.ones(n[1]+1)) + self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)] + elif(self.dim == 3): + area1 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], vh[2]))) + area2 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + area3 = np.outer(vh[0], Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)] + return self._area - def area(): - doc = "Construct face areas of the 3D model as 1d array." - - def fget(self): - if(self._area is None): - # Ensure that we are working with column vectors - vh = self.h - # The number of cell centers in each direction - n = self.vnC - # Compute areas of cell faces - if(self.dim == 1): - self._area = np.ones(n[0]+1) - elif(self.dim == 2): - area1 = np.outer(np.ones(n[0]+1), vh[1]) - area2 = np.outer(vh[0], np.ones(n[1]+1)) - self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)] - elif(self.dim == 3): - area1 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], vh[2]))) - area2 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) - area3 = np.outer(vh[0], Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) - self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)] - return self._area - return locals() - _area = None - area = property(**area()) - - def edge(): - doc = "Construct edge legnths of the 3D model as 1d array." - - def fget(self): - if(self._edge is None): - # Ensure that we are working with column vectors - vh = self.h - # The number of cell centers in each direction - n = self.vnC - # Compute edge lengths - if(self.dim == 1): - self._edge = Utils.mkvc(vh[0]) - elif(self.dim == 2): - l1 = np.outer(vh[0], np.ones(n[1]+1)) - l2 = np.outer(np.ones(n[0]+1), vh[1]) - self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2)] - elif(self.dim == 3): - l1 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) - l2 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) - l3 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) - self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2), Utils.mkvc(l3)] - return self._edge - return locals() - _edge = None - edge = property(**edge()) + @property + def edge(self): + """Construct edge legnths of the 3D model as 1d array.""" + if getattr(self, '_edge', None) is None: + # Ensure that we are working with column vectors + vh = self.h + # The number of cell centers in each direction + n = self.vnC + # Compute edge lengths + if(self.dim == 1): + self._edge = Utils.mkvc(vh[0]) + elif(self.dim == 2): + l1 = np.outer(vh[0], np.ones(n[1]+1)) + l2 = np.outer(np.ones(n[0]+1), vh[1]) + self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2)] + elif(self.dim == 3): + l1 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) + l2 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + l3 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2), Utils.mkvc(l3)] + return self._edge # --------------- Methods --------------------- diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py index 47b30243..b92cf0eb 100644 --- a/SimPEG/Mesh/__init__.py +++ b/SimPEG/Mesh/__init__.py @@ -1,8 +1,8 @@ -from Cyl1DMesh import Cyl1DMesh -from TensorMesh import TensorMesh -from TreeMesh import TreeMesh -from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh from BaseMesh import BaseMesh, BaseRectangularMesh +from TensorMesh import TensorMesh +from CylMesh import CylMesh +from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh +from TreeMesh import TreeMesh from TensorView import TensorView from LomView import LomView from InnerProducts import InnerProducts diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py new file mode 100644 index 00000000..5ce86e0c --- /dev/null +++ b/SimPEG/Tests/test_cylMesh.py @@ -0,0 +1,29 @@ +import unittest +import sys +from SimPEG import * + + +class TestCylMesh(unittest.TestCase): + + def setUp(self): + hx = np.ones(10) + hz = np.ones(5) + + self.mesh = Mesh.CylMesh([hx, 1,hz]) + + def test_cylMeshInheritance(self): + self.assertTrue(isinstance(self.mesh, Mesh.BaseMesh)) + + def test_cylMeshDimensions(self): + self.assertTrue(self.mesh.dim == 3) + + def test_cylMesh_nc(self): + self.assertTrue(np.all(self.mesh.vnC == [10, 1, 5])) + + def test_cylMesh_nc_xyz(self): + self.assertTrue(self.mesh.nCx == 10) + self.assertTrue(self.mesh.nCy == 1) + self.assertTrue(self.mesh.nCz == 5) + +if __name__ == '__main__': + unittest.main()