From 250108277f9c1914f6e47acc4f3008bd7218543d Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 14 Feb 2014 17:33:27 -0800 Subject: [PATCH 01/23] Initial clean up of the CylMesh code. Inherits from TensorMesh. Many more things to do! --- SimPEG/Mesh/BaseMesh.py | 14 +- SimPEG/Mesh/Cyl1DMesh.py | 476 ----------------------------------- SimPEG/Mesh/CylMesh.py | 337 +++++++++++++++++++++++++ SimPEG/Mesh/TensorMesh.py | 344 +++++++++++-------------- SimPEG/Mesh/__init__.py | 8 +- SimPEG/Tests/test_cylMesh.py | 29 +++ 6 files changed, 520 insertions(+), 688 deletions(-) delete mode 100644 SimPEG/Mesh/Cyl1DMesh.py create mode 100644 SimPEG/Mesh/CylMesh.py create mode 100644 SimPEG/Tests/test_cylMesh.py 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() From dc563ff6c5ce932042a9c7a8aab7cb328aa09d4e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 14 Feb 2014 20:25:57 -0800 Subject: [PATCH 02/23] Correct counting on cylinder mesh! --- SimPEG/Mesh/CylMesh.py | 70 +++++++++++++++++++++++++++--------- SimPEG/Tests/test_cylMesh.py | 43 ++++++++++++++++------ 2 files changed, 87 insertions(+), 26 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 594c739c..3f01be59 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -31,6 +31,15 @@ class CylMesh(TensorMesh): TensorMesh.__init__(self, h, x0) + @property + def nNx(self): + """ + Number of nodes in the x-direction + + :rtype: int + :return: nNx + """ + return self.nCx @property def nNy(self): @@ -40,7 +49,7 @@ class CylMesh(TensorMesh): :rtype: int :return: nNy """ - return self.nCy + return self.nCy - 1 @property def nN(self): @@ -49,15 +58,29 @@ class CylMesh(TensorMesh): :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 nFx(self): + """ + Number of x-faces + + :rtype: int + :return: nFx + """ + return self.nC + + @property + def vnFx(self): + """ + Number of x-faces in each direction + + :rtype: numpy.array (dim, ) + :return: vnFx + """ + return self.vnC + @property def nFy(self): """ @@ -69,19 +92,34 @@ class CylMesh(TensorMesh): 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]) + def nEx(self): + """ + Number of x-edges + + :rtype: int + :return: nEx + """ + return (self._n + np.r_[0,-1,1]).prod() @property - def nF(self): - """Total number of faces""" - return self.nFr + self.vnFz + def nEy(self): + """ + Number of y-edges + + :rtype: int + :return: nEy + """ + return (self._n + np.r_[0,0,1]).prod() @property - def nE(self): - """Number of edges""" - return self.nN + def nEz(self): + """ + Number of z-edges + + :rtype: int + :return: nEz + """ + return (self._n + np.r_[0,-1,0]).prod() @property def vectorNx(self): diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 5ce86e0c..aa1af471 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -3,12 +3,11 @@ import sys from SimPEG import * -class TestCylMesh(unittest.TestCase): +class TestCyl1DMesh(unittest.TestCase): def setUp(self): - hx = np.ones(10) - hz = np.ones(5) - + hx = np.ones(3) + hz = np.ones(2) self.mesh = Mesh.CylMesh([hx, 1,hz]) def test_cylMeshInheritance(self): @@ -17,13 +16,37 @@ class TestCylMesh(unittest.TestCase): 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) + def test_cylMesh_numbers(self): + self.assertTrue(self.mesh.nCx == 3) self.assertTrue(self.mesh.nCy == 1) - self.assertTrue(self.mesh.nCz == 5) + self.assertTrue(self.mesh.nCz == 2) + self.assertTrue(np.all(self.mesh.vnC == [3, 1, 2])) + + self.assertTrue(self.mesh.nN == 0) + self.assertTrue(self.mesh.nNx == 3) + self.assertTrue(self.mesh.nNy == 0) + self.assertTrue(self.mesh.nNz == 3) + self.assertTrue(np.all(self.mesh.vnN == [3, 0, 3])) + + self.assertTrue(self.mesh.nFx == 6) + self.assertTrue(np.all(self.mesh.vnFx == [3, 1, 2])) + self.assertTrue(self.mesh.nFy == 0) + self.assertTrue(np.all(self.mesh.vnFy == [3, 0, 2])) + self.assertTrue(self.mesh.nFz == 9) + self.assertTrue(np.all(self.mesh.vnFz == [3, 1, 3])) + self.assertTrue(self.mesh.nF == 15) + self.assertTrue(np.all(self.mesh.vnF == [6, 0, 9])) + + self.assertTrue(self.mesh.nEx == 0) + self.assertTrue(np.all(self.mesh.vnEx == [3, 0, 3])) + self.assertTrue(self.mesh.nEy == 9) + self.assertTrue(np.all(self.mesh.vnEy == [3, 1, 3])) + self.assertTrue(self.mesh.nEz == 0) + self.assertTrue(np.all(self.mesh.vnEz == [3, 0, 2])) + self.assertTrue(self.mesh.nE == 9) + + + if __name__ == '__main__': unittest.main() From 1129fc1c66d3779c1e7d400774a78247801fd8f9 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 14 Feb 2014 21:44:46 -0800 Subject: [PATCH 03/23] vectors in CylMesh --- SimPEG/Mesh/CylMesh.py | 41 +++++++++++++++++++---- SimPEG/Tests/test_cylMesh.py | 65 ++++++++++++++++++++++++++++++++---- 2 files changed, 93 insertions(+), 13 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 3f01be59..935c8ce0 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -6,7 +6,7 @@ from TensorMesh import TensorMesh class CylMesh(TensorMesh): """ - CylMesh is a mesh class for cylindrically problems + CylMesh is a mesh class for cylindrical problems """ _meshType = 'CYL' @@ -15,7 +15,8 @@ class CylMesh(TensorMesh): 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" + assert type(x0) == np.ndarray, "x0 must be an ndarray" + assert x0.size == 3, "x0 must have 3 elements" else: x0 = np.r_[0, 0, 0] @@ -39,7 +40,9 @@ class CylMesh(TensorMesh): :rtype: int :return: nNx """ - return self.nCx + if self.nCy == 1: + return self.nCx + return self.nCx + 1 @property def nNy(self): @@ -49,7 +52,9 @@ class CylMesh(TensorMesh): :rtype: int :return: nNy """ - return self.nCy - 1 + if self.nCy == 1: + return self.nCy - 1 + return self.nCy @property def nN(self): @@ -121,16 +126,38 @@ class CylMesh(TensorMesh): """ return (self._n + np.r_[0,-1,0]).prod() + @property + def vectorCCx(self): + """Cell-centered grid vector (1D) in the x direction.""" + if self.nCy == 1: + return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 - self.hx[0]/2 + return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + + @property + def vectorCCy(self): + """Cell-centered grid vector (1D) in the y direction.""" + return np.r_[0, self.hy[:-1]] + @property def vectorNx(self): - """Nodal grid vector (1D) in the r direction""" - return self.hr.cumsum() + """Nodal grid vector (1D) in the x direction.""" + if self.nCy == 1: + return self.hx.cumsum() - self.hx[0]/2 + return np.r_[0, self.hx].cumsum() + + @property + def vectorNy(self): + """Nodal grid vector (1D) in the y direction.""" + return np.r_[0, self.hy[:-1].cumsum()] + self.hy[0]*0.5 @property def edge(self): """Edge lengths""" if getattr(self, '_edge', None) is None: - self._edge = 2*pi*self.gridN[:,0] + if self.nCy == 1: + self._edge = 2*pi*self.gridN[:,0] + else: + raise NotImplementedError('edges not implemented for 3D cyl mesh') return self._edge @property diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index aa1af471..c0db7982 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -3,20 +3,18 @@ import sys from SimPEG import * -class TestCyl1DMesh(unittest.TestCase): +class TestCyl2DMesh(unittest.TestCase): def setUp(self): - hx = np.ones(3) - hz = np.ones(2) + hx = np.r_[1,1,0.5] + hz = np.r_[2,1] 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_numbers(self): + self.assertTrue(self.mesh.dim == 3) self.assertTrue(self.mesh.nCx == 3) self.assertTrue(self.mesh.nCy == 1) self.assertTrue(self.mesh.nCz == 2) @@ -45,7 +43,62 @@ class TestCyl1DMesh(unittest.TestCase): self.assertTrue(np.all(self.mesh.vnEz == [3, 0, 2])) self.assertTrue(self.mesh.nE == 9) + def test_vectorsCC(self): + v = np.r_[0, 1, 1.75] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0) + v = np.r_[0] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCy)) == 0) + v = np.r_[1, 2.5] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCz)) == 0) + def test_vectorsN(self): + v = np.r_[0.5, 1.5, 2] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNx)) == 0) + v = np.r_[np.pi] #This is kinda a fake. But it is where it would be if there was a radial connection + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNy)) == 0) + v = np.r_[0, 2, 3.] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNz)) == 0) + + def test_dimensions(self): + v = np.r_[0.5, 1.5, 2, 0.5, 1.5, 2, 0.5, 1.5, 2] * 2 * np.pi + self.assertTrue(np.linalg.norm((v-self.mesh.edge)) == 0) + + +class TestCyl3DMesh(unittest.TestCase): + + def setUp(self): + hx = np.r_[1,1,0.5] + hy = np.r_[np.pi, np.pi] + hz = np.r_[2,1] + self.mesh = Mesh.CylMesh([hx, hy,hz]) + + def test_cylMesh_numbers(self): + self.assertTrue(self.mesh.nCx == 3) + self.assertTrue(self.mesh.nCy == 2) + self.assertTrue(self.mesh.nCz == 2) + self.assertTrue(np.all(self.mesh.vnC == [3, 2, 2])) + + self.assertTrue(self.mesh.nN == 24) + self.assertTrue(self.mesh.nNx == 4) + self.assertTrue(self.mesh.nNy == 2) + self.assertTrue(self.mesh.nNz == 3) + self.assertTrue(np.all(self.mesh.vnN == [4, 2, 3])) + + def test_vectorsCC(self): + v = np.r_[0.5, 1.5, 2.25] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0) + v = np.r_[0, np.pi] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCy)) == 0) + v = np.r_[1, 2.5] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCz)) == 0) + + def test_vectorsN(self): + v = np.r_[0, 1, 2, 2.5] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNx)) == 0) + v = np.r_[np.pi/2, 1.5*np.pi] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNy)) == 0) + v = np.r_[0, 2, 3] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNz)) == 0) if __name__ == '__main__': From 9bea954dff876273fe3c5049ff697478e7423048 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 15 Feb 2014 14:08:41 -0800 Subject: [PATCH 04/23] vol and area tests passing --- SimPEG/Mesh/CylMesh.py | 14 +++++++++----- SimPEG/Mesh/InnerProducts.py | 2 ++ SimPEG/Tests/test_cylMesh.py | 18 ++++++++++++++++-- 3 files changed, 27 insertions(+), 7 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 935c8ce0..34f19c1e 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -157,15 +157,17 @@ class CylMesh(TensorMesh): if self.nCy == 1: self._edge = 2*pi*self.gridN[:,0] else: - raise NotImplementedError('edges not implemented for 3D cyl mesh') + raise NotImplementedError('edges not yet implemented for 3D cyl mesh') 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)) + if self.nCy > 1: + raise NotImplementedError('area not yet implemented for 3D cyl mesh') + areaR = np.kron(self.hz, 2*pi*self.vectorNx) + areaZ = np.kron(np.ones_like(self.vectorNz),pi*(self.vectorNx**2 - np.r_[0, self.vectorNx[:-1]]**2)) self._area = np.r_[areaR, areaZ] return self._area @@ -173,7 +175,9 @@ class CylMesh(TensorMesh): 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) + if self.nCy > 1: + raise NotImplementedError('vol not yet implemented for 3D cyl mesh') + az = pi*(self.vectorNx**2 - np.r_[0, self.vectorNx[:-1]]**2) self._vol = np.kron(self.hz,az) return self._vol @@ -295,7 +299,7 @@ class CylMesh(TensorMesh): loc = np.atleast_2d(loc) - assert np.all(loc[:,0]<=self.vectorNr.max()) & \ + assert np.all(loc[:,0]<=self.vectorNx.max()) & \ np.all(loc[:,1]>=self.vectorNz.min()) & \ np.all(loc[:,1]<=self.vectorNz.max()), \ "Points outside of mesh" diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 427e5ee2..a5257c8d 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -303,6 +303,8 @@ class InnerProducts(object): def _makeTensor(M, sigma): if sigma is None: # default is ones sigma = np.ones((M.nC, 1)) + elif type(sigma) is float: + sigma = np.ones(self.nC)*sigma if M.dim == 2: if sigma.size == M.nC: # Isotropic! diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index c0db7982..3103e1d5 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -15,6 +15,8 @@ class TestCyl2DMesh(unittest.TestCase): def test_cylMesh_numbers(self): self.assertTrue(self.mesh.dim == 3) + + self.assertTrue(self.mesh.nC == 6) self.assertTrue(self.mesh.nCx == 3) self.assertTrue(self.mesh.nCy == 1) self.assertTrue(self.mesh.nCz == 2) @@ -42,6 +44,7 @@ class TestCyl2DMesh(unittest.TestCase): self.assertTrue(self.mesh.nEz == 0) self.assertTrue(np.all(self.mesh.vnEz == [3, 0, 2])) self.assertTrue(self.mesh.nE == 9) + self.assertTrue(np.all(self.mesh.vnE == [0, 9, 0])) def test_vectorsCC(self): v = np.r_[0, 1, 1.75] @@ -60,9 +63,20 @@ class TestCyl2DMesh(unittest.TestCase): self.assertTrue(np.linalg.norm((v-self.mesh.vectorNz)) == 0) def test_dimensions(self): - v = np.r_[0.5, 1.5, 2, 0.5, 1.5, 2, 0.5, 1.5, 2] * 2 * np.pi - self.assertTrue(np.linalg.norm((v-self.mesh.edge)) == 0) + edge = np.r_[0.5, 1.5, 2, 0.5, 1.5, 2, 0.5, 1.5, 2] * 2 * np.pi + self.assertTrue(np.linalg.norm((edge-self.mesh.edge)) == 0) + r = np.r_[0, 0.5, 1.5, 2] + a = r[1:]*2*np.pi + areaX = np.r_[2*a,a] + a = (r[1:]**2 - r[:-1]**2)*np.pi + areaZ = np.r_[a,a,a] + area = np.r_[areaX, areaZ] + self.assertTrue(np.linalg.norm((area-self.mesh.area)) == 0) + + a = (r[1:]**2 - r[:-1]**2)*np.pi + vol = np.r_[2*a,a] + self.assertTrue(np.linalg.norm((vol-self.mesh.vol)) == 0) class TestCyl3DMesh(unittest.TestCase): From 2bdc3958b5a54bb3d3dd525b5d937cf6a8c1d120 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 16 Feb 2014 14:42:43 -0800 Subject: [PATCH 05/23] Split out base tensor mesh stuff into different mixing class. --- SimPEG/Mesh/CylMesh.py | 57 +++++--- SimPEG/Mesh/TensorMesh.py | 252 +++++++++++++++++---------------- SimPEG/Tests/test_operators.py | 1 + 3 files changed, 169 insertions(+), 141 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 34f19c1e..eb557735 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -2,35 +2,21 @@ 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 +from TensorMesh import BaseTensorMesh -class CylMesh(TensorMesh): +class CylMesh(BaseTensorMesh): """ CylMesh is a mesh class for cylindrical problems """ _meshType = 'CYL' + _unitDimensions = [1, 2*np.pi, 1] + 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 type(x0) == np.ndarray, "x0 must be an ndarray" - assert x0.size == 3, "x0 must have 3 elements" - 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) + 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" @property def nNx(self): @@ -185,6 +171,35 @@ class CylMesh(TensorMesh): # Operators #################################################### + @property + def faceDiv(self): + """Construct divergence operator (face-stg to cell-centres).""" + raise NotImplementedError('faceDiv not yet implemented') + @property + def faceDivx(self): + """Construct divergence operator in the x component (face-stg to cell-centres).""" + raise NotImplementedError('faceDivx not yet implemented') + + @property + def faceDivy(self): + """Construct divergence operator in the y component (face-stg to cell-centres).""" + raise NotImplementedError('faceDivy not yet implemented') + + @property + def faceDivz(self): + """Construct divergence operator in the z component (face-stg to cell-centres).""" + raise NotImplementedError('faceDivz not yet implemented') + + @property + def nodalGrad(self): + """Construct gradient operator (nodes to edges).""" + raise NotImplementedError('nodalGrad not yet implemented') + + @property + def nodalLaplacian(self): + """Construct laplacian operator (nodes to edges).""" + raise NotImplementedError('nodalLaplacian not yet implemented') + @property def edgeCurl(self): """The edgeCurl property.""" diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 693359bf..6f292d0f 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -4,38 +4,13 @@ from TensorView import TensorView from DiffOperators import DiffOperators from InnerProducts import InnerProducts -class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): - """ - TensorMesh is a mesh class that deals with tensor product meshes. - - Any Mesh that has a constant width along the entire axis - such that it can defined by a single width vector, called 'h'. - - :: - - hx = np.array([1,1,1]) - hy = np.array([1,2]) - hz = np.array([1,1,1,1]) - - mesh = Mesh.TensorMesh([hx, hy, hz]) - - Example of a padded tensor mesh: - - .. plot:: - - from SimPEG import Mesh, Utils - M = Mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) - M.plotGrid() - - For a quick tensor mesh on a (10x12x15) unit cube:: - - mesh = Mesh.TensorMesh([10, 12, 15]) - - """ +class BaseTensorMesh(BaseRectangularMesh): __metaclass__ = Utils.SimPEGMetaClass - _meshType = 'TENSOR' + _meshType = 'BASETENSOR' + + _unitDimensions = [1, 1, 1] def __init__(self, h_in, x0=None): assert type(h_in) is list, 'h_in must be a list' @@ -43,7 +18,7 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): for i, h_i in enumerate(h_in): if type(h_i) in [int, long, float]: # This gives you something over the unit cube. - h_i = np.ones(int(h_i))/int(h_i) + h_i = self._unitDimensions[i] * 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. @@ -54,57 +29,6 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): # Ensure h contains 1D vectors self._h = [Utils.mkvc(x.astype(float)) for x in h] - def __str__(self): - outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) - def printH(hx, outStr=''): - i = -1 - while True: - i = i + 1 - if i > hx.size: - break - elif i == hx.size: - break - h = hx[i] - n = 1 - for j in range(i+1, hx.size): - if hx[j] == h: - n = n + 1 - i = i + 1 - else: - break - - if n == 1: - outStr = outStr + ' {0:.2f},'.format(h) - else: - outStr = outStr + ' {0:d}*{1:.2f},'.format(n,h) - - return outStr[:-1] - - if self.dim == 1: - outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) - outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) - outStr = outStr + printH(self.hx, outStr='\n hx:') - pass - elif self.dim == 2: - outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) - outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1]) - outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) - outStr = outStr + '\n nCy: {0:d}'.format(self.nCy) - outStr = outStr + printH(self.hx, outStr='\n hx:') - outStr = outStr + printH(self.hy, outStr='\n hy:') - elif self.dim == 3: - outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) - outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1]) - outStr = outStr + '\n z0: {0:.2f}'.format(self.x0[2]) - outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) - outStr = outStr + '\n nCy: {0:d}'.format(self.nCy) - outStr = outStr + '\n nCz: {0:d}'.format(self.nCz) - outStr = outStr + printH(self.hx, outStr='\n hx:') - outStr = outStr + printH(self.hy, outStr='\n hy:') - outStr = outStr + printH(self.hz, outStr='\n hz:') - - return outStr - @property def h(self): """h is a list containing the cell widths of the tensor mesh in each dimension.""" @@ -211,6 +135,133 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): self._gridEz = Utils.ndgrid(self.getTensor('Ez')) return self._gridEz + def getTensor(self, locType): + """ Returns a tensor list. + + :param str locType: What tensor (see below) + :rtype: list + :return: list of the tensors that make up the mesh. + + locType can be:: + + 'Ex' -> x-component of field defined on edges + 'Ey' -> y-component of field defined on edges + 'Ez' -> z-component of field defined on edges + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces + 'N' -> scalar field defined on nodes + 'CC' -> scalar field defined on cell centers + """ + + if locType is 'Fx': + ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] + elif locType is 'Fy': + ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] + elif locType is 'Fz': + ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] + elif locType is 'Ex': + ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] + elif locType is 'Ey': + ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] + elif locType is 'Ez': + ten = [self.vectorNx , self.vectorNy , self.vectorCCz] + elif locType is 'CC': + ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] + elif locType is 'N': + ten = [self.vectorNx , self.vectorNy , self.vectorNz ] + + return [t for t in ten if t is not None] + + +class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): + """ + TensorMesh is a mesh class that deals with tensor product meshes. + + Any Mesh that has a constant width along the entire axis + such that it can defined by a single width vector, called 'h'. + + :: + + hx = np.array([1,1,1]) + hy = np.array([1,2]) + hz = np.array([1,1,1,1]) + + mesh = Mesh.TensorMesh([hx, hy, hz]) + + Example of a padded tensor mesh: + + .. plot:: + + from SimPEG import Mesh, Utils + M = Mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + M.plotGrid() + + For a quick tensor mesh on a (10x12x15) unit cube:: + + mesh = Mesh.TensorMesh([10, 12, 15]) + + """ + + __metaclass__ = Utils.SimPEGMetaClass + + _meshType = 'TENSOR' + + def __init__(self, h_in, x0=None): + BaseTensorMesh.__init__(self, h_in, x0) + + def __str__(self): + outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) + def printH(hx, outStr=''): + i = -1 + while True: + i = i + 1 + if i > hx.size: + break + elif i == hx.size: + break + h = hx[i] + n = 1 + for j in range(i+1, hx.size): + if hx[j] == h: + n = n + 1 + i = i + 1 + else: + break + + if n == 1: + outStr = outStr + ' {0:.2f},'.format(h) + else: + outStr = outStr + ' {0:d}*{1:.2f},'.format(n,h) + + return outStr[:-1] + + if self.dim == 1: + outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) + outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) + outStr = outStr + printH(self.hx, outStr='\n hx:') + pass + elif self.dim == 2: + outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) + outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1]) + outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) + outStr = outStr + '\n nCy: {0:d}'.format(self.nCy) + outStr = outStr + printH(self.hx, outStr='\n hx:') + outStr = outStr + printH(self.hy, outStr='\n hy:') + elif self.dim == 3: + outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) + outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1]) + outStr = outStr + '\n z0: {0:.2f}'.format(self.x0[2]) + outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) + outStr = outStr + '\n nCy: {0:d}'.format(self.nCy) + outStr = outStr + '\n nCz: {0:d}'.format(self.nCz) + outStr = outStr + printH(self.hx, outStr='\n hx:') + outStr = outStr + printH(self.hy, outStr='\n hy:') + outStr = outStr + printH(self.hz, outStr='\n hz:') + + return outStr + + # --------------- Geometries --------------------- @property def vol(self): @@ -274,45 +325,6 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): # --------------- Methods --------------------- - def getTensor(self, locType): - """ Returns a tensor list. - - :param str locType: What tensor (see below) - :rtype: list - :return: list of the tensors that make up the mesh. - - locType can be:: - - 'Ex' -> x-component of field defined on edges - 'Ey' -> y-component of field defined on edges - 'Ez' -> z-component of field defined on edges - 'Fx' -> x-component of field defined on faces - 'Fy' -> y-component of field defined on faces - 'Fz' -> z-component of field defined on faces - 'N' -> scalar field defined on nodes - 'CC' -> scalar field defined on cell centers - """ - - if locType is 'Fx': - ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] - elif locType is 'Fy': - ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] - elif locType is 'Fz': - ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] - elif locType is 'Ex': - ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] - elif locType is 'Ey': - ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] - elif locType is 'Ez': - ten = [self.vectorNx , self.vectorNy , self.vectorCCz] - elif locType is 'CC': - ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] - elif locType is 'N': - ten = [self.vectorNx , self.vectorNy , self.vectorNz ] - - return [t for t in ten if t is not None] - - def isInside(self, pts): """ Determines if a set of points are inside a mesh. diff --git a/SimPEG/Tests/test_operators.py b/SimPEG/Tests/test_operators.py index 67741aa7..76884b96 100644 --- a/SimPEG/Tests/test_operators.py +++ b/SimPEG/Tests/test_operators.py @@ -3,6 +3,7 @@ import unittest from TestUtils import OrderTest import matplotlib.pyplot as plt +#TODO: 'randomTensorMesh' MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM'] call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) From 7dd8f9e83e594db6cdceb05d22a4fa2322c3b2f1 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 16 Feb 2014 15:32:20 -0800 Subject: [PATCH 06/23] faceDiv (untested) --- SimPEG/Mesh/CylMesh.py | 39 +++++++++++++++++++++++++++++++----- SimPEG/Mesh/DiffOperators.py | 6 +++--- SimPEG/Tests/test_cylMesh.py | 3 +++ 3 files changed, 40 insertions(+), 8 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index eb557735..9ab89e64 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -1,7 +1,7 @@ import numpy as np import scipy.sparse as sp from scipy.constants import pi -from SimPEG.Utils import mkvc, ndgrid, sdiag +from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, ddx, av, avExtrap from TensorMesh import BaseTensorMesh class CylMesh(BaseTensorMesh): @@ -174,21 +174,50 @@ class CylMesh(BaseTensorMesh): @property def faceDiv(self): """Construct divergence operator (face-stg to cell-centres).""" - raise NotImplementedError('faceDiv not yet implemented') + if getattr(self, '_faceDiv', None) is None: + n = self.vnC + # Compute faceDivergence operator on faces + D1 = self.faceDivx + D3 = self.faceDivz + if self.nCy == 1: + D = sp.hstack((D1, D3), format="csr") + elif self.nCy > 1: + D2 = self.faceDivy + D = sp.hstack((D1, D2, D3), format="csr") + self._faceDiv = D + return self._faceDiv + @property def faceDivx(self): """Construct divergence operator in the x component (face-stg to cell-centres).""" - raise NotImplementedError('faceDivx not yet implemented') + if getattr(self, '_faceDivx', None) is None: + D1 = kron3(speye(self.nCz), speye(self.nCy), ddx(self.nCx)[:,1:]) + S = self.r(self.area, 'F', 'Fx', 'V') + V = self.vol + self._faceDivx = sdiag(1/V)*D1*sdiag(S) + return self._faceDivx @property def faceDivy(self): """Construct divergence operator in the y component (face-stg to cell-centres).""" - raise NotImplementedError('faceDivy not yet implemented') + raise NotImplementedError('Wrapping the ddx is not yet implemented.') + if getattr(self, '_faceDivy', None) is None: + # TODO: this needs to wrap to join up faces which are connected in the cylinder + D2 = kron3(speye(self.nCz), ddx(self.nCy), speye(self.nCx)) + S = self.r(self.area, 'F', 'Fy', 'V') + V = self.vol + self._faceDivy = sdiag(1/V)*D2*sdiag(S) + return self._faceDivy @property def faceDivz(self): """Construct divergence operator in the z component (face-stg to cell-centres).""" - raise NotImplementedError('faceDivz not yet implemented') + if getattr(self, '_faceDivz', None) is None: + D3 = kron3(ddx(self.nCz), speye(self.nCy), speye(self.nCx)) + S = self.r(self.area, 'F', 'Fz', 'V') + V = self.vol + self._faceDivz = sdiag(1/V)*D3*sdiag(S) + return self._faceDivz @property def nodalGrad(self): diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 1712964e..f23171c9 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -167,7 +167,7 @@ class DiffOperators(object): elif(self.dim == 3): D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0])) # Compute areas of cell faces & volumes - S = self.r(self.area, 'F','Fx', 'V') + S = self.r(self.area, 'F', 'Fx', 'V') V = self.vol self._faceDivx = sdiag(1/V)*D1*sdiag(S) @@ -190,7 +190,7 @@ class DiffOperators(object): elif(self.dim == 3): D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0])) # Compute areas of cell faces & volumes - S = self.r(self.area, 'F','Fy', 'V') + S = self.r(self.area, 'F', 'Fy', 'V') V = self.vol self._faceDivy = sdiag(1/V)*D2*sdiag(S) @@ -210,7 +210,7 @@ class DiffOperators(object): # Compute faceDivergence operator on faces D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0])) # Compute areas of cell faces & volumes - S = self.r(self.area, 'F','Fz', 'V') + S = self.r(self.area, 'F', 'Fz', 'V') V = self.vol self._faceDivz = sdiag(1/V)*D3*sdiag(S) diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 3103e1d5..b0adf9bf 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -78,6 +78,9 @@ class TestCyl2DMesh(unittest.TestCase): vol = np.r_[2*a,a] self.assertTrue(np.linalg.norm((vol-self.mesh.vol)) == 0) + def test_faceDiv(self): + print self.mesh.faceDiv + class TestCyl3DMesh(unittest.TestCase): def setUp(self): From d5c09ce03686911271b87120d36b8f3a0d1134b2 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 16 Feb 2014 15:43:49 -0800 Subject: [PATCH 07/23] Combine view functions into one file. --- SimPEG/Mesh/LogicallyOrthogonalMesh.py | 2 +- SimPEG/Mesh/LomView.py | 104 ------------------------- SimPEG/Mesh/TensorMesh.py | 2 +- SimPEG/Mesh/{TensorView.py => View.py} | 98 +++++++++++++++++++++++ SimPEG/Mesh/__init__.py | 5 -- SimPEG/Tests/test_basemesh.py | 2 +- SimPEG/Tests/test_cylMesh.py | 3 - 7 files changed, 101 insertions(+), 115 deletions(-) delete mode 100644 SimPEG/Mesh/LomView.py rename SimPEG/Mesh/{TensorView.py => View.py} (80%) diff --git a/SimPEG/Mesh/LogicallyOrthogonalMesh.py b/SimPEG/Mesh/LogicallyOrthogonalMesh.py index 480c4f06..49196e23 100644 --- a/SimPEG/Mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/Mesh/LogicallyOrthogonalMesh.py @@ -2,7 +2,7 @@ from SimPEG import Utils, np from BaseMesh import BaseRectangularMesh from DiffOperators import DiffOperators from InnerProducts import InnerProducts -from LomView import LomView +from View import LomView # Some helper functions. length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5 diff --git a/SimPEG/Mesh/LomView.py b/SimPEG/Mesh/LomView.py deleted file mode 100644 index 4eceb918..00000000 --- a/SimPEG/Mesh/LomView.py +++ /dev/null @@ -1,104 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import matplotlib -from mpl_toolkits.mplot3d import Axes3D -from SimPEG.Utils import mkvc - - -class LomView(object): - """ - Provides viewing functions for LogicallyOrthogonalMesh - - This class is inherited by LogicallyOrthogonalMesh - - """ - def __init__(self): - pass - - def plotGrid(self, length=0.05, showIt=False): - """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. - - - .. plot:: - :include-source: - - from SimPEG import Mesh, Utils - X, Y = Utils.exampleLomGird([3,3],'rotate') - M = Mesh.LogicallyOrthogonalMesh([X, Y]) - M.plotGrid(showIt=True) - - """ - NN = self.r(self.gridN, 'N', 'N', 'M') - if self.dim == 2: - fig = plt.figure(2) - fig.clf() - ax = plt.subplot(111) - X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten() - - X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten() - - X = np.r_[X1, X2] - Y = np.r_[Y1, Y2] - - plt.plot(X, Y) - - plt.hold(True) - Nx = self.r(self.normals, 'F', 'Fx', 'V') - Ny = self.r(self.normals, 'F', 'Fy', 'V') - Tx = self.r(self.tangents, 'E', 'Ex', 'V') - Ty = self.r(self.tangents, 'E', 'Ey', 'V') - - plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo') - - nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten() - nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten() - plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs') - plt.plot(nX, nY, 'r-') - - nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten() - nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten() - #plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs') - plt.plot(nX, nY, 'g-') - - tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten() - tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten() - plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^') - plt.plot(tX, tY, 'r-') - - nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten() - nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten() - #plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^') - plt.plot(nX, nY, 'g-') - plt.axis('equal') - - elif self.dim == 3: - fig = plt.figure(3) - fig.clf() - ax = fig.add_subplot(111, projection='3d') - X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten() - Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten() - - X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten() - Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten() - - X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten() - Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten() - Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten() - - X = np.r_[X1, X2, X3] - Y = np.r_[Y1, Y2, Y3] - Z = np.r_[Z1, Z2, Z3] - - plt.plot(X, Y, 'b', zs=Z) - ax.set_zlabel('x3') - - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - ax.set_ylabel('x2') - - if showIt: plt.show() diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 6f292d0f..9717affa 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -1,6 +1,6 @@ from SimPEG import Utils, np, sp from BaseMesh import BaseRectangularMesh -from TensorView import TensorView +from View import TensorView from DiffOperators import DiffOperators from InnerProducts import InnerProducts diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/View.py similarity index 80% rename from SimPEG/Mesh/TensorView.py rename to SimPEG/Mesh/View.py index e4dd9243..c6b1d6c5 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/View.py @@ -426,3 +426,101 @@ class TensorView(object): return animate(fig, animateFrame, frames=len(frames)) + +class LomView(object): + """ + Provides viewing functions for LogicallyOrthogonalMesh + + This class is inherited by LogicallyOrthogonalMesh + + """ + def __init__(self): + pass + + def plotGrid(self, length=0.05, showIt=False): + """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. + + + .. plot:: + :include-source: + + from SimPEG import Mesh, Utils + X, Y = Utils.exampleLomGird([3,3],'rotate') + M = Mesh.LogicallyOrthogonalMesh([X, Y]) + M.plotGrid(showIt=True) + + """ + NN = self.r(self.gridN, 'N', 'N', 'M') + if self.dim == 2: + fig = plt.figure(2) + fig.clf() + ax = plt.subplot(111) + X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten() + Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten() + + X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten() + Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten() + + X = np.r_[X1, X2] + Y = np.r_[Y1, Y2] + + plt.plot(X, Y) + + plt.hold(True) + Nx = self.r(self.normals, 'F', 'Fx', 'V') + Ny = self.r(self.normals, 'F', 'Fy', 'V') + Tx = self.r(self.tangents, 'E', 'Ex', 'V') + Ty = self.r(self.tangents, 'E', 'Ey', 'V') + + plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo') + + nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten() + nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten() + plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs') + plt.plot(nX, nY, 'r-') + + nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten() + nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten() + #plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs') + plt.plot(nX, nY, 'g-') + + tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten() + tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten() + plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^') + plt.plot(tX, tY, 'r-') + + nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten() + nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten() + #plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^') + plt.plot(nX, nY, 'g-') + plt.axis('equal') + + elif self.dim == 3: + fig = plt.figure(3) + fig.clf() + ax = fig.add_subplot(111, projection='3d') + X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten() + Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten() + Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten() + + X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten() + Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten() + Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten() + + X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten() + Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten() + Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten() + + X = np.r_[X1, X2, X3] + Y = np.r_[Y1, Y2, Y3] + Z = np.r_[Z1, Z2, Z3] + + plt.plot(X, Y, 'b', zs=Z) + ax.set_zlabel('x3') + + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + + if showIt: plt.show() diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py index b92cf0eb..de46f675 100644 --- a/SimPEG/Mesh/__init__.py +++ b/SimPEG/Mesh/__init__.py @@ -1,9 +1,4 @@ -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 -from DiffOperators import DiffOperators diff --git a/SimPEG/Tests/test_basemesh.py b/SimPEG/Tests/test_basemesh.py index eb3acc21..d482505e 100644 --- a/SimPEG/Tests/test_basemesh.py +++ b/SimPEG/Tests/test_basemesh.py @@ -1,6 +1,6 @@ import unittest import sys -from SimPEG.Mesh import BaseRectangularMesh +from SimPEG.Mesh.BaseMesh import BaseRectangularMesh import numpy as np diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index b0adf9bf..b6006bff 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -10,9 +10,6 @@ class TestCyl2DMesh(unittest.TestCase): hz = np.r_[2,1] self.mesh = Mesh.CylMesh([hx, 1,hz]) - def test_cylMeshInheritance(self): - self.assertTrue(isinstance(self.mesh, Mesh.BaseMesh)) - def test_cylMesh_numbers(self): self.assertTrue(self.mesh.dim == 3) From 8b5ca46ecf0f9aab6bb16fc0473bf9819b6f9220 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 16 Feb 2014 18:05:11 -0800 Subject: [PATCH 08/23] started testing of faceDiv --- SimPEG/Mesh/CylMesh.py | 46 +++++++++++++++++++++++++++++ SimPEG/Tests/TestUtils.py | 21 +++++++++---- SimPEG/Tests/test_cylMesh.py | 57 ++++++++++++++++++++++++++++++++++-- 3 files changed, 117 insertions(+), 7 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 9ab89e64..da0927e5 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -136,6 +136,45 @@ class CylMesh(BaseTensorMesh): """Nodal grid vector (1D) in the y direction.""" return np.r_[0, self.hy[:-1].cumsum()] + self.hy[0]*0.5 + + def getTensor(self, locType): + """ Returns a tensor list. + + :param str locType: What tensor (see below) + :rtype: list + :return: list of the tensors that make up the mesh. + + locType can be:: + + 'Ex' -> x-component of field defined on edges + 'Ey' -> y-component of field defined on edges + 'Ez' -> z-component of field defined on edges + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces + 'N' -> scalar field defined on nodes + 'CC' -> scalar field defined on cell centers + """ + + if locType is 'Fx': + ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] + elif locType is 'Fy': + ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] + elif locType is 'Fz': + ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] + elif locType is 'Ex': + ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] + elif locType is 'Ey': + ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] + elif locType is 'Ez': + ten = [self.vectorNx , self.vectorNy , self.vectorCCz] + elif locType is 'CC': + ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] + elif locType is 'N': + ten = [self.vectorNx , self.vectorNy , self.vectorNz ] + + return [t for t in ten if t is not None] + @property def edge(self): """Edge lengths""" @@ -219,6 +258,13 @@ class CylMesh(BaseTensorMesh): self._faceDivz = sdiag(1/V)*D3*sdiag(S) return self._faceDivz + + @property + def cellGrad(self): + """The cell centered Gradient, takes you to cell faces.""" + raise NotImplementedError('Cell Grad is not yet implemented.') + + @property def nodalGrad(self): """Construct gradient operator (nodes to edges).""" diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index d84408d4..a980546e 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt from numpy.linalg import norm from SimPEG.Utils import mkvc, sdiag from SimPEG import Utils -from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh +from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh, CylMesh import numpy as np import scipy.sparse as sp import unittest @@ -88,10 +88,7 @@ class OrderTest(unittest.TestCase): """ if 'TensorMesh' in self._meshType: if 'uniform' in self._meshType: - h1 = np.ones(nc)/nc - h2 = np.ones(nc)/nc - h3 = np.ones(nc)/nc - h = [h1, h2, h3] + h = [nc, nc, nc] elif 'random' in self._meshType: h1 = np.random.rand(nc) h2 = np.random.rand(nc) @@ -104,6 +101,20 @@ class OrderTest(unittest.TestCase): max_h = max([np.max(hi) for hi in self.M.h]) return max_h + elif 'CylMesh' in self._meshType: + if 'uniform' in self._meshType: + h = [nc, nc, nc] + else: + raise Exception('Unexpected meshType') + + if self.meshDimension == 2: + self.M = CylMesh([h[0], 1, h[2]]) + max_h = max([np.max(hi) for hi in [self.M.hx, self.M.hz]]) + elif self.meshDimension == 3: + self.M = CylMesh(h) + max_h = max([np.max(hi) for hi in self.M.h]) + return max_h + elif 'LOM' in self._meshType: if 'uniform' in self._meshType: kwrd = 'rect' diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index b6006bff..584bcc2d 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -1,6 +1,7 @@ import unittest import sys from SimPEG import * +from TestUtils import OrderTest class TestCyl2DMesh(unittest.TestCase): @@ -75,8 +76,60 @@ class TestCyl2DMesh(unittest.TestCase): vol = np.r_[2*a,a] self.assertTrue(np.linalg.norm((vol-self.mesh.vol)) == 0) - def test_faceDiv(self): - print self.mesh.faceDiv + def test_gridSizes(self): + self.assertTrue(self.mesh.gridCC.shape == (self.mesh.nC, 3)) + # self.assertTrue(self.mesh.gridN.shape == (self.mesh.nN, 3)) + + self.assertTrue(self.mesh.gridFx.shape == (self.mesh.nFx, 3)) + # self.assertTrue(self.mesh.gridFy.shape == (self.mesh.nFy, 3)) + self.assertTrue(self.mesh.gridFz.shape == (self.mesh.nFz, 3)) + + # self.assertTrue(self.mesh.gridEx.shape == (self.mesh.nEx, 3)) + self.assertTrue(self.mesh.gridEy.shape == (self.mesh.nEy, 3)) + # self.assertTrue(self.mesh.gridEz.shape == (self.mesh.nEz, 3)) + + +MESHTYPES = ['uniformCylMesh'] +MESHDIMENSION = 2 +call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) +call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) +cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] +cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)] +cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy))) +cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey))) +cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz))) +cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez))) + + +class TestFaceDiv(OrderTest): + name = "FaceDiv" + meshTypes = MESHTYPES + meshDimension = MESHDIMENSION + + def getError(self): + + funX = lambda x, y, z: np.cos(2*np.pi*y) + funY = lambda x, y, z: 1+x*0 + funZ = lambda x, y, z: np.cos(2*np.pi*x) + + solX = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*z) + solY = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x) + solZ = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*y) + + Ec = cartE3(self.M, funX, funY, funZ) + print Ec.shape, self.M.nE + print self.M + E = self.M.projectEdgeVector(Ec) + + Fc = cartF3(self.M, solX, solY, solZ) + curlE_anal = self.M.projectFaceVector(Fc) + + curlE = self.M.edgeCurl.dot(E) + err = np.linalg.norm((curlE - curlE_anal), np.inf) + return err + + # def test_order(self): + # self.orderTest() class TestCyl3DMesh(unittest.TestCase): From ba52e16bea6759d9ad4a0003a13071d64ae07662 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 11:04:00 -0800 Subject: [PATCH 09/23] updated where cell centers are and how I count the first h in 2D case --- SimPEG/Mesh/CylMesh.py | 6 ++++-- SimPEG/Tests/test_cylMesh.py | 10 +++++----- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index da0927e5..376cb789 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -116,7 +116,7 @@ class CylMesh(BaseTensorMesh): def vectorCCx(self): """Cell-centered grid vector (1D) in the x direction.""" if self.nCy == 1: - return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 - self.hx[0]/2 + return np.r_[-self.hx[0]*0.5, self.hx[:-1].cumsum()] + self.hx*0.5 # - self.hx[0]/2 return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 @property @@ -128,12 +128,14 @@ class CylMesh(BaseTensorMesh): def vectorNx(self): """Nodal grid vector (1D) in the x direction.""" if self.nCy == 1: - return self.hx.cumsum() - self.hx[0]/2 + return self.hx.cumsum()# - self.hx[0]/2 return np.r_[0, self.hx].cumsum() @property def vectorNy(self): """Nodal grid vector (1D) in the y direction.""" + if self.nCy == 1: + return np.r_[0] return np.r_[0, self.hy[:-1].cumsum()] + self.hy[0]*0.5 diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 584bcc2d..ebde658c 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -45,7 +45,7 @@ class TestCyl2DMesh(unittest.TestCase): self.assertTrue(np.all(self.mesh.vnE == [0, 9, 0])) def test_vectorsCC(self): - v = np.r_[0, 1, 1.75] + v = np.r_[0, 1.5, 2.25] self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0) v = np.r_[0] self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCy)) == 0) @@ -53,18 +53,18 @@ class TestCyl2DMesh(unittest.TestCase): self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCz)) == 0) def test_vectorsN(self): - v = np.r_[0.5, 1.5, 2] + v = np.r_[1, 2, 2.5] self.assertTrue(np.linalg.norm((v-self.mesh.vectorNx)) == 0) - v = np.r_[np.pi] #This is kinda a fake. But it is where it would be if there was a radial connection + v = np.r_[0] self.assertTrue(np.linalg.norm((v-self.mesh.vectorNy)) == 0) v = np.r_[0, 2, 3.] self.assertTrue(np.linalg.norm((v-self.mesh.vectorNz)) == 0) def test_dimensions(self): - edge = np.r_[0.5, 1.5, 2, 0.5, 1.5, 2, 0.5, 1.5, 2] * 2 * np.pi + edge = np.r_[1, 2, 2.5, 1, 2, 2.5, 1, 2, 2.5] * 2 * np.pi self.assertTrue(np.linalg.norm((edge-self.mesh.edge)) == 0) - r = np.r_[0, 0.5, 1.5, 2] + r = np.r_[0, 1, 2, 2.5] a = r[1:]*2*np.pi areaX = np.r_[2*a,a] a = (r[1:]**2 - r[:-1]**2)*np.pi From 0fb229a3c0a43da8955708cca339ee86f0f47c26 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 11:26:14 -0800 Subject: [PATCH 10/23] update tensors and tests --- SimPEG/Mesh/CylMesh.py | 50 +++++------------------------------- SimPEG/Mesh/TensorMesh.py | 6 +++++ SimPEG/Tests/test_cylMesh.py | 37 +++++++++++++++++++++++--- 3 files changed, 46 insertions(+), 47 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 376cb789..226b8720 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -115,9 +115,8 @@ class CylMesh(BaseTensorMesh): @property def vectorCCx(self): """Cell-centered grid vector (1D) in the x direction.""" - if self.nCy == 1: - return np.r_[-self.hx[0]*0.5, self.hx[:-1].cumsum()] + self.hx*0.5 # - self.hx[0]/2 - return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + firstEl = -self.hx[0]*0.5 if self.nCy == 1 else 0 + return np.r_[firstEl, self.hx[:-1].cumsum()] + self.hx*0.5 @property def vectorCCy(self): @@ -128,55 +127,18 @@ class CylMesh(BaseTensorMesh): def vectorNx(self): """Nodal grid vector (1D) in the x direction.""" if self.nCy == 1: - return self.hx.cumsum()# - self.hx[0]/2 + return self.hx.cumsum() return np.r_[0, self.hx].cumsum() @property def vectorNy(self): """Nodal grid vector (1D) in the y direction.""" if self.nCy == 1: + # There aren't really any nodes, but all the grids need + # somewhere to live, why not zero?! return np.r_[0] return np.r_[0, self.hy[:-1].cumsum()] + self.hy[0]*0.5 - - def getTensor(self, locType): - """ Returns a tensor list. - - :param str locType: What tensor (see below) - :rtype: list - :return: list of the tensors that make up the mesh. - - locType can be:: - - 'Ex' -> x-component of field defined on edges - 'Ey' -> y-component of field defined on edges - 'Ez' -> z-component of field defined on edges - 'Fx' -> x-component of field defined on faces - 'Fy' -> y-component of field defined on faces - 'Fz' -> z-component of field defined on faces - 'N' -> scalar field defined on nodes - 'CC' -> scalar field defined on cell centers - """ - - if locType is 'Fx': - ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] - elif locType is 'Fy': - ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] - elif locType is 'Fz': - ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] - elif locType is 'Ex': - ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] - elif locType is 'Ey': - ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] - elif locType is 'Ez': - ten = [self.vectorNx , self.vectorNy , self.vectorCCz] - elif locType is 'CC': - ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] - elif locType is 'N': - ten = [self.vectorNx , self.vectorNy , self.vectorNz ] - - return [t for t in ten if t is not None] - @property def edge(self): """Edge lengths""" @@ -270,6 +232,8 @@ class CylMesh(BaseTensorMesh): @property def nodalGrad(self): """Construct gradient operator (nodes to edges).""" + if self.nCy == 1: + raise Exception('Nodal grad does not make sense for cylindrically symmetric mesh.') raise NotImplementedError('nodalGrad not yet implemented') @property diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 9a68d12a..44623883 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -95,6 +95,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridFx(self): + if self.nFx == 0: return """Face staggered grid in the x direction.""" if getattr(self, '_gridFx', None) is None: self._gridFx = Utils.ndgrid(self.getTensor('Fx')) @@ -102,6 +103,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridFy(self): + if self.nFy == 0: return """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')) @@ -109,6 +111,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridFz(self): + if self.nFz == 0: return """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')) @@ -116,6 +119,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridEx(self): + if self.nEx == 0: return """Edge staggered grid in the x direction.""" if getattr(self, '_gridEx', None) is None: self._gridEx = Utils.ndgrid(self.getTensor('Ex')) @@ -123,6 +127,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridEy(self): + if self.nEy == 0: return """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')) @@ -130,6 +135,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridEz(self): + if self.nEz == 0: return """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')) diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index ebde658c..9eaf0998 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -78,15 +78,44 @@ class TestCyl2DMesh(unittest.TestCase): def test_gridSizes(self): self.assertTrue(self.mesh.gridCC.shape == (self.mesh.nC, 3)) - # self.assertTrue(self.mesh.gridN.shape == (self.mesh.nN, 3)) + self.assertTrue(self.mesh.gridN.shape == (9, 3)) self.assertTrue(self.mesh.gridFx.shape == (self.mesh.nFx, 3)) - # self.assertTrue(self.mesh.gridFy.shape == (self.mesh.nFy, 3)) + self.assertTrue(self.mesh.gridFy is None) self.assertTrue(self.mesh.gridFz.shape == (self.mesh.nFz, 3)) - # self.assertTrue(self.mesh.gridEx.shape == (self.mesh.nEx, 3)) + self.assertTrue(self.mesh.gridEx is None) self.assertTrue(self.mesh.gridEy.shape == (self.mesh.nEy, 3)) - # self.assertTrue(self.mesh.gridEz.shape == (self.mesh.nEz, 3)) + self.assertTrue(self.mesh.gridEz is None) + + def test_gridCC(self): + x = np.r_[0,1.5,2.25,0,1.5,2.25] + y = np.zeros(6) + z = np.r_[1,1,1,2.5,2.5,2.5] + G = np.c_[x,y,z] + self.assertTrue(np.linalg.norm((G-self.mesh.gridCC).ravel()) == 0) + + def test_gridFx(self): + x = np.r_[1,2,2.5,1,2,2.5] + y = np.zeros(6) + z = np.r_[1,1,1,2.5,2.5,2.5] + G = np.c_[x,y,z] + self.assertTrue(np.linalg.norm((G-self.mesh.gridFx).ravel()) == 0) + + def test_gridFz(self): + x = np.r_[0,1.5,2.25,0,1.5,2.25,0,1.5,2.25] + y = np.zeros(9) + z = np.r_[0,0,0,2,2,2,3,3,3.] + G = np.c_[x,y,z] + self.assertTrue(np.linalg.norm((G-self.mesh.gridFz).ravel()) == 0) + + def test_gridEy(self): + x = np.r_[1,2,2.5,1,2,2.5,1,2,2.5] + y = np.zeros(9) + z = np.r_[0,0,0,2,2,2,3,3,3.] + G = np.c_[x,y,z] + self.assertTrue(np.linalg.norm((G-self.mesh.gridEy).ravel()) == 0) + MESHTYPES = ['uniformCylMesh'] From c859735be3eec11310104bc6f08dce17c0dfe2ca Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 11:27:15 -0800 Subject: [PATCH 11/23] test gridN --- SimPEG/Tests/test_cylMesh.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 9eaf0998..86e47475 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -95,6 +95,13 @@ class TestCyl2DMesh(unittest.TestCase): G = np.c_[x,y,z] self.assertTrue(np.linalg.norm((G-self.mesh.gridCC).ravel()) == 0) + def test_gridN(self): + x = np.r_[1,2,2.5,1,2,2.5,1,2,2.5] + y = np.zeros(9) + z = np.r_[0,0,0,2,2,2,3,3,3.] + G = np.c_[x,y,z] + self.assertTrue(np.linalg.norm((G-self.mesh.gridN).ravel()) == 0) + def test_gridFx(self): x = np.r_[1,2,2.5,1,2,2.5] y = np.zeros(6) From 9f7cf7eed10730965aac94ab2910966e194f8333 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 11:28:51 -0800 Subject: [PATCH 12/23] test edge area and volume --- SimPEG/Tests/test_cylMesh.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 86e47475..70c34c88 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -60,10 +60,11 @@ class TestCyl2DMesh(unittest.TestCase): v = np.r_[0, 2, 3.] self.assertTrue(np.linalg.norm((v-self.mesh.vectorNz)) == 0) - def test_dimensions(self): + def test_edge(self): edge = np.r_[1, 2, 2.5, 1, 2, 2.5, 1, 2, 2.5] * 2 * np.pi self.assertTrue(np.linalg.norm((edge-self.mesh.edge)) == 0) + def test_area(self): r = np.r_[0, 1, 2, 2.5] a = r[1:]*2*np.pi areaX = np.r_[2*a,a] @@ -72,6 +73,8 @@ class TestCyl2DMesh(unittest.TestCase): area = np.r_[areaX, areaZ] self.assertTrue(np.linalg.norm((area-self.mesh.area)) == 0) + def test_vol(self): + r = np.r_[0, 1, 2, 2.5] a = (r[1:]**2 - r[:-1]**2)*np.pi vol = np.r_[2*a,a] self.assertTrue(np.linalg.norm((vol-self.mesh.vol)) == 0) From 8f8cbb1c19c7152a53359b9510eb74df7040af53 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 12:08:08 -0800 Subject: [PATCH 13/23] Move some counting down to base mesh. --- SimPEG/Mesh/BaseMesh.py | 88 ++++++++++++++++++++++++++++++++++++ SimPEG/Mesh/CylMesh.py | 71 ++++------------------------- SimPEG/Tests/test_cylMesh.py | 35 +++++++++++++- 3 files changed, 131 insertions(+), 63 deletions(-) diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py index cecf6d7b..9af638b8 100644 --- a/SimPEG/Mesh/BaseMesh.py +++ b/SimPEG/Mesh/BaseMesh.py @@ -399,6 +399,94 @@ class BaseRectangularMesh(BaseMesh): """ return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None]) + ################################## + # Redo the numbering so they are dependent of the vector numbers + ################################## + + @property + def nC(self): + """ + Total number of cells + + :rtype: int + :return: nC + """ + return self.vnC.prod() + + @property + def nN(self): + """ + Total number of nodes + + :rtype: int + :return: nN + """ + return self.vnN.prod() + + @property + def nEx(self): + """ + Number of x-edges + + :rtype: int + :return: nEx + """ + return self.vnEx.prod() + + @property + def nEy(self): + """ + Number of y-edges + + :rtype: int + :return: nEy + """ + if self.dim < 2: return + return self.vnEy.prod() + + @property + def nEz(self): + """ + Number of z-edges + + :rtype: int + :return: nEz + """ + if self.dim < 3: return + return self.vnEz.prod() + + @property + def nFx(self): + """ + Number of x-faces + + :rtype: int + :return: nFx + """ + return self.vnFx.prod() + + @property + def nFy(self): + """ + Number of y-faces + + :rtype: int + :return: nFy + """ + if self.dim < 2: return + return self.vnFy.prod() + + @property + def nFz(self): + """ + Number of z-faces + + :rtype: int + :return: nFz + """ + if self.dim < 3: return + return self.vnFz.prod() + def r(self, x, xType='CC', outType='CC', format='V'): """ Mesh.r is a quick reshape command that will do the best it can at giving you what you want. diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 226b8720..0502ed31 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -26,8 +26,7 @@ class CylMesh(BaseTensorMesh): :rtype: int :return: nNx """ - if self.nCy == 1: - return self.nCx + if self.nCy == 1: return self.nCx return self.nCx + 1 @property @@ -38,30 +37,9 @@ class CylMesh(BaseTensorMesh): :rtype: int :return: nNy """ - if self.nCy == 1: - return self.nCy - 1 + if self.nCy == 1: return 0 return self.nCy - @property - def nN(self): - """ - Total number of nodes - - :rtype: int - :return: nN - """ - return (np.r_[self.nNx, self.nNy, self.nNz]).prod() - - @property - def nFx(self): - """ - Number of x-faces - - :rtype: int - :return: nFx - """ - return self.nC - @property def vnFx(self): """ @@ -73,44 +51,15 @@ class CylMesh(BaseTensorMesh): return self.vnC @property - def nFy(self): + def vnEy(self): """ - Number of y-faces + Number of y-edges in each direction - :rtype: int - :return: nFy + :rtype: numpy.array (dim, ) + :return: vnEy or None if dim < 2 """ - return (self.vnC + np.r_[0,-1,0][:self.dim]).prod() - - @property - def nEx(self): - """ - Number of x-edges - - :rtype: int - :return: nEx - """ - return (self._n + np.r_[0,-1,1]).prod() - - @property - def nEy(self): - """ - Number of y-edges - - :rtype: int - :return: nEy - """ - return (self._n + np.r_[0,0,1]).prod() - - @property - def nEz(self): - """ - Number of z-edges - - :rtype: int - :return: nEz - """ - return (self._n + np.r_[0,-1,0]).prod() + nNx = self.nNx if self.nCy == 1 else self.nNx - 1 + return np.r_[nNx, self.nCy, self.nNz] @property def vectorCCx(self): @@ -232,8 +181,8 @@ class CylMesh(BaseTensorMesh): @property def nodalGrad(self): """Construct gradient operator (nodes to edges).""" - if self.nCy == 1: - raise Exception('Nodal grad does not make sense for cylindrically symmetric mesh.') + # Nodal grad does not make sense for cylindrically symmetric mesh. + if self.nCy == 1: return None raise NotImplementedError('nodalGrad not yet implemented') @property diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 70c34c88..18db1686 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -11,21 +11,24 @@ class TestCyl2DMesh(unittest.TestCase): hz = np.r_[2,1] self.mesh = Mesh.CylMesh([hx, 1,hz]) - def test_cylMesh_numbers(self): + def test_dim(self): self.assertTrue(self.mesh.dim == 3) + def test_nC(self): self.assertTrue(self.mesh.nC == 6) self.assertTrue(self.mesh.nCx == 3) self.assertTrue(self.mesh.nCy == 1) self.assertTrue(self.mesh.nCz == 2) self.assertTrue(np.all(self.mesh.vnC == [3, 1, 2])) + def test_nN(self): self.assertTrue(self.mesh.nN == 0) self.assertTrue(self.mesh.nNx == 3) self.assertTrue(self.mesh.nNy == 0) self.assertTrue(self.mesh.nNz == 3) self.assertTrue(np.all(self.mesh.vnN == [3, 0, 3])) + def test_nF(self): self.assertTrue(self.mesh.nFx == 6) self.assertTrue(np.all(self.mesh.vnFx == [3, 1, 2])) self.assertTrue(self.mesh.nFy == 0) @@ -35,6 +38,7 @@ class TestCyl2DMesh(unittest.TestCase): self.assertTrue(self.mesh.nF == 15) self.assertTrue(np.all(self.mesh.vnF == [6, 0, 9])) + def test_nE(self): self.assertTrue(self.mesh.nEx == 0) self.assertTrue(np.all(self.mesh.vnEx == [3, 0, 3])) self.assertTrue(self.mesh.nEy == 9) @@ -126,6 +130,9 @@ class TestCyl2DMesh(unittest.TestCase): G = np.c_[x,y,z] self.assertTrue(np.linalg.norm((G-self.mesh.gridEy).ravel()) == 0) + def test_lightOperators(self): + self.assertTrue(self.mesh.nodalGrad is None) + MESHTYPES = ['uniformCylMesh'] @@ -178,18 +185,42 @@ class TestCyl3DMesh(unittest.TestCase): hz = np.r_[2,1] self.mesh = Mesh.CylMesh([hx, hy,hz]) - def test_cylMesh_numbers(self): + def test_dim(self): + self.assertTrue(self.mesh.dim == 3) + + def test_nC(self): self.assertTrue(self.mesh.nCx == 3) self.assertTrue(self.mesh.nCy == 2) self.assertTrue(self.mesh.nCz == 2) self.assertTrue(np.all(self.mesh.vnC == [3, 2, 2])) + def test_nN(self): self.assertTrue(self.mesh.nN == 24) self.assertTrue(self.mesh.nNx == 4) self.assertTrue(self.mesh.nNy == 2) self.assertTrue(self.mesh.nNz == 3) self.assertTrue(np.all(self.mesh.vnN == [4, 2, 3])) + def test_nF(self): + self.assertTrue(self.mesh.nFx == 12) + self.assertTrue(np.all(self.mesh.vnFx == [3, 2, 2])) + self.assertTrue(self.mesh.nFy == 12) + self.assertTrue(np.all(self.mesh.vnFy == [3, 2, 2])) + self.assertTrue(self.mesh.nFz == 18) + self.assertTrue(np.all(self.mesh.vnFz == [3, 2, 3])) + self.assertTrue(self.mesh.nF == 42) + self.assertTrue(np.all(self.mesh.vnF == [12, 12, 18])) + + def test_nE(self): + self.assertTrue(self.mesh.nEx == 18) + self.assertTrue(np.all(self.mesh.vnEx == [3, 2, 3])) + self.assertTrue(self.mesh.nEy == 18) + self.assertTrue(np.all(self.mesh.vnEy == [3, 2, 3])) + # self.assertTrue(self.mesh.nEz == 0) + # self.assertTrue(np.all(self.mesh.vnEz == [3, 0, 2])) + # self.assertTrue(self.mesh.nE == 9) + # self.assertTrue(np.all(self.mesh.vnE == [0, 9, 0])) + def test_vectorsCC(self): v = np.r_[0.5, 1.5, 2.25] self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0) From 4d7d10dcd3da6d127f230b3f38c9e10c804ed295 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 13:42:15 -0800 Subject: [PATCH 14/23] updates to number of edges in the z direction --- SimPEG/Mesh/CylMesh.py | 27 ++++++++++++++++++++++++++- SimPEG/Tests/test_cylMesh.py | 8 ++++---- 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 0502ed31..482d5724 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -61,6 +61,31 @@ class CylMesh(BaseTensorMesh): nNx = self.nNx if self.nCy == 1 else self.nNx - 1 return np.r_[nNx, self.nCy, self.nNz] + @property + def vnEz(self): + """ + Number of z-edges in each direction + + :rtype: numpy.array (dim, ) + :return: vnEz or None if nCy > 1 + """ + if self.nCy == 1: + return np.r_[self.nNx, self.nNy, self.nCz] + else: + return None + + @property + def nEz(self): + """ + Number of z-edges + + :rtype: int + :return: nEz + """ + if self.nCy == 1: + return self.vnEz.prod() + return (np.r_[self.nNx-1, self.nNy, self.nCz]).prod() + self.nCz + @property def vectorCCx(self): """Cell-centered grid vector (1D) in the x direction.""" @@ -116,7 +141,7 @@ class CylMesh(BaseTensorMesh): if self.nCy > 1: raise NotImplementedError('vol not yet implemented for 3D cyl mesh') az = pi*(self.vectorNx**2 - np.r_[0, self.vectorNx[:-1]]**2) - self._vol = np.kron(self.hz,az) + self._vol = np.kron(self.hz, az) return self._vol #################################################### diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 18db1686..1be073ef 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -216,10 +216,10 @@ class TestCyl3DMesh(unittest.TestCase): self.assertTrue(np.all(self.mesh.vnEx == [3, 2, 3])) self.assertTrue(self.mesh.nEy == 18) self.assertTrue(np.all(self.mesh.vnEy == [3, 2, 3])) - # self.assertTrue(self.mesh.nEz == 0) - # self.assertTrue(np.all(self.mesh.vnEz == [3, 0, 2])) - # self.assertTrue(self.mesh.nE == 9) - # self.assertTrue(np.all(self.mesh.vnE == [0, 9, 0])) + self.assertTrue(self.mesh.nEz == 12 + 2) + self.assertTrue(self.mesh.vnEz is None) + self.assertTrue(self.mesh.nE == 50) + self.assertTrue(np.all(self.mesh.vnE == [18, 18, 14])) def test_vectorsCC(self): v = np.r_[0.5, 1.5, 2.25] From 897b43bf80c4da17989c53fb1d66ef3cabd43307 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 14:05:37 -0800 Subject: [PATCH 15/23] refactor tensor grids --- SimPEG/Mesh/TensorMesh.py | 50 ++++++++++++++++----------------------- 1 file changed, 20 insertions(+), 30 deletions(-) diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 44623883..2b50103e 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -4,6 +4,12 @@ from View import TensorView from DiffOperators import DiffOperators from InnerProducts import InnerProducts +def _getTensorGrid(self, key): + if getattr(self, '_grid' + key, None) is None: + setattr(self, '_grid' + key, Utils.ndgrid(self.getTensor(key))) + return getattr(self, '_grid' + key) + + class BaseTensorMesh(BaseRectangularMesh): __metaclass__ = Utils.SimPEGMetaClass @@ -82,64 +88,48 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridCC(self): """Cell-centered grid.""" - if getattr(self, '_gridCC', None) is None: - self._gridCC = Utils.ndgrid(self.getTensor('CC')) - return self._gridCC + return _getTensorGrid(self, 'CC') @property def gridN(self): """Nodal grid.""" - if getattr(self, '_gridN', None) is None: - self._gridN = Utils.ndgrid(self.getTensor('N')) - return self._gridN + return _getTensorGrid(self, 'N') @property def gridFx(self): - if self.nFx == 0: return """Face staggered grid in the x direction.""" - if getattr(self, '_gridFx', None) is None: - self._gridFx = Utils.ndgrid(self.getTensor('Fx')) - return self._gridFx + if self.nFx == 0: return + return _getTensorGrid(self, 'Fx') @property def gridFy(self): - if self.nFy == 0: return """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 + if self.nFy == 0 or self.dim < 2: return + return _getTensorGrid(self, 'Fy') @property def gridFz(self): - if self.nFz == 0: return """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 + if self.nFz == 0 or self.dim < 3: return + return _getTensorGrid(self, 'Fz') @property def gridEx(self): - if self.nEx == 0: return """Edge staggered grid in the x direction.""" - if getattr(self, '_gridEx', None) is None: - self._gridEx = Utils.ndgrid(self.getTensor('Ex')) - return self._gridEx + if self.nEx == 0: return + return _getTensorGrid(self, 'Ex') @property def gridEy(self): - if self.nEy == 0: return """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 + if self.nEy == 0 or self.dim < 2: return + return _getTensorGrid(self, 'Ey') @property def gridEz(self): - if self.nEz == 0: return """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 + if self.nEz == 0 or self.dim < 3: return + return _getTensorGrid(self, 'Ez') def getTensor(self, locType): """ Returns a tensor list. From fee512dda393b24456d254d8fa3d2d912333ea16 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 15:31:31 -0800 Subject: [PATCH 16/23] changes to CC at middle of first h_r --- SimPEG/Mesh/CylMesh.py | 3 +-- SimPEG/Mesh/TensorMesh.py | 2 +- SimPEG/Tests/test_cylMesh.py | 48 +++++++++++++++++------------------- 3 files changed, 24 insertions(+), 29 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 482d5724..21a178c5 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -89,8 +89,7 @@ class CylMesh(BaseTensorMesh): @property def vectorCCx(self): """Cell-centered grid vector (1D) in the x direction.""" - firstEl = -self.hx[0]*0.5 if self.nCy == 1 else 0 - return np.r_[firstEl, self.hx[:-1].cumsum()] + self.hx*0.5 + return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 @property def vectorCCy(self): diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 2b50103e..8204d5ae 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -22,7 +22,7 @@ class BaseTensorMesh(BaseRectangularMesh): assert type(h_in) is list, 'h_in must be a list' h = range(len(h_in)) for i, h_i in enumerate(h_in): - if type(h_i) in [int, long, float]: + if type(h_i) in [int, long, float, np.int_]: # This gives you something over the unit cube. h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i) assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i) diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 1be073ef..bbe44bb7 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -49,7 +49,7 @@ class TestCyl2DMesh(unittest.TestCase): self.assertTrue(np.all(self.mesh.vnE == [0, 9, 0])) def test_vectorsCC(self): - v = np.r_[0, 1.5, 2.25] + v = np.r_[0.5, 1.5, 2.25] self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0) v = np.r_[0] self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCy)) == 0) @@ -96,7 +96,7 @@ class TestCyl2DMesh(unittest.TestCase): self.assertTrue(self.mesh.gridEz is None) def test_gridCC(self): - x = np.r_[0,1.5,2.25,0,1.5,2.25] + x = np.r_[0.5,1.5,2.25,0.5,1.5,2.25] y = np.zeros(6) z = np.r_[1,1,1,2.5,2.5,2.5] G = np.c_[x,y,z] @@ -117,7 +117,7 @@ class TestCyl2DMesh(unittest.TestCase): self.assertTrue(np.linalg.norm((G-self.mesh.gridFx).ravel()) == 0) def test_gridFz(self): - x = np.r_[0,1.5,2.25,0,1.5,2.25,0,1.5,2.25] + x = np.r_[0.5,1.5,2.25,0.5,1.5,2.25,0.5,1.5,2.25] y = np.zeros(9) z = np.r_[0,0,0,2,2,2,3,3,3.] G = np.c_[x,y,z] @@ -137,14 +137,14 @@ class TestCyl2DMesh(unittest.TestCase): MESHTYPES = ['uniformCylMesh'] MESHDIMENSION = 2 -call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) +call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 2]) call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) -cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] -cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)] -cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy))) -cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey))) -cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz))) -cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez))) +cyl_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] +cyl_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)] +cylF2 = lambda M, fx, fy: np.vstack((cyl_row2(M.gridFx, fx, fy), cyl_row2(M.gridFz, fx, fy))) +# cylE2 = lambda M, ex, ey: np.vstack((cyl_row2(M.gridEx, ex, ey), cyl_row2(M.gridEy, ex, ey))) +# cylF3 = lambda M, fx, fy, fz: np.vstack((cyl_row3(M.gridFx, fx, fy, fz), cyl_row3(M.gridFy, fx, fy, fz), cyl_row3(M.gridFz, fx, fy, fz))) +# cylE3 = lambda M, ex, ey, ez: np.vstack((cyl_row3(M.gridEx, ex, ey, ez), cyl_row3(M.gridEy, ex, ey, ez), cyl_row3(M.gridEz, ex, ey, ez))) class TestFaceDiv(OrderTest): @@ -154,28 +154,24 @@ class TestFaceDiv(OrderTest): def getError(self): - funX = lambda x, y, z: np.cos(2*np.pi*y) - funY = lambda x, y, z: 1+x*0 - funZ = lambda x, y, z: np.cos(2*np.pi*x) + funR = lambda r, z: np.sin(2.*np.pi*r) + funZ = lambda r, z: np.sin(2.*np.pi*z) - solX = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*z) - solY = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x) - solZ = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*y) + sol = lambda r, t, z: (2*np.pi*r*np.cos(2*np.pi*r) + np.sin(2*np.pi*r))/r + 2*np.pi*np.cos(2*np.pi*z) - Ec = cartE3(self.M, funX, funY, funZ) - print Ec.shape, self.M.nE - print self.M - E = self.M.projectEdgeVector(Ec) + Fc = cylF2(self.M, funR, funZ) + Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]] + F = self.M.projectFaceVector(Fc) - Fc = cartF3(self.M, solX, solY, solZ) - curlE_anal = self.M.projectFaceVector(Fc) + divF = self.M.faceDiv.dot(F) + tm = Mesh.TensorMesh([self.M.nCx, self.M.nCz]) + divF_anal = call3(sol, self.M.gridCC) - curlE = self.M.edgeCurl.dot(E) - err = np.linalg.norm((curlE - curlE_anal), np.inf) + err = np.linalg.norm((divF-divF_anal), np.inf) return err - # def test_order(self): - # self.orderTest() + def test_order(self): + self.orderTest() class TestCyl3DMesh(unittest.TestCase): From 72510420a6dff1de407b157e3046619718c603df Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 23 Feb 2014 13:52:23 -0800 Subject: [PATCH 17/23] testCylEdgeCurl start --- SimPEG/Mesh/CylMesh.py | 2 ++ SimPEG/Tests/test_cylMesh.py | 17 ++++++++++++++--- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 21a178c5..267533a6 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -217,6 +217,8 @@ class CylMesh(BaseTensorMesh): @property def edgeCurl(self): """The edgeCurl property.""" + if self.nCy > 1: + raise NotImplementedError('Edge curl not yet implemented for nCy > 1') 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") diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index bbe44bb7..cd6ba6da 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -136,7 +136,6 @@ class TestCyl2DMesh(unittest.TestCase): MESHTYPES = ['uniformCylMesh'] -MESHDIMENSION = 2 call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 2]) call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) cyl_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] @@ -147,10 +146,10 @@ cylF2 = lambda M, fx, fy: np.vstack((cyl_row2(M.gridFx, fx, fy), cyl_row2(M.grid # cylE3 = lambda M, ex, ey, ez: np.vstack((cyl_row3(M.gridEx, ex, ey, ez), cyl_row3(M.gridEy, ex, ey, ez), cyl_row3(M.gridEz, ex, ey, ez))) -class TestFaceDiv(OrderTest): +class TestFaceDiv2D(OrderTest): name = "FaceDiv" meshTypes = MESHTYPES - meshDimension = MESHDIMENSION + meshDimension = 2 def getError(self): @@ -173,6 +172,18 @@ class TestFaceDiv(OrderTest): def test_order(self): self.orderTest() +class TestEdgeCurl2D(OrderTest): + name = "EdgeCurl" + meshTypes = MESHTYPES + meshDimension = 2 + + def getError(self): + + #TODO! + + # def test_order(self): + # self.orderTest() + class TestCyl3DMesh(unittest.TestCase): def setUp(self): From 5bcdc5d46acb32f90b3454b7399d9195233c703d Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 6 Mar 2014 19:05:58 -0800 Subject: [PATCH 18/23] rearranging of interpolation code. --- SimPEG/Mesh/CylMesh.py | 207 ++++++++++++++++++----------------- SimPEG/Mesh/TensorMesh.py | 159 ++++++++++++++------------- SimPEG/Tests/test_cylMesh.py | 1 - 3 files changed, 186 insertions(+), 181 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 267533a6..c44afebe 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -4,6 +4,37 @@ from scipy.constants import pi from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, ddx, av, avExtrap from TensorMesh import BaseTensorMesh + + # 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 + class CylMesh(BaseTensorMesh): """ CylMesh is a mesh class for cylindrical problems @@ -313,125 +344,97 @@ class CylMesh(BaseTensorMesh): """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 + # 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 + # :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:: + # 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 - """ + # '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) + # loc = np.atleast_2d(loc) - assert np.all(loc[:,0]<=self.vectorNx.max()) & \ - np.all(loc[:,1]>=self.vectorNz.min()) & \ - np.all(loc[:,1]<=self.vectorNz.max()), \ - "Points outside of mesh" + # assert np.all(loc[:,0]<=self.vectorNx.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) + # 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) + # 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 + # 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,:] + # 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 + # 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,:] + # 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] + # 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] + # 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,:] + # 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 + # 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 + # 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() + # 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 377c948b..cc3bcbb1 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -170,6 +170,87 @@ class BaseTensorMesh(BaseRectangularMesh): return [t for t in ten if t is not None] + # --------------- Methods --------------------- + + def isInside(self, pts, locType='N'): + """ + Determines if a set of points are inside a mesh. + + :param numpy.ndarray pts: Location of points to test + :rtype numpy.ndarray + :return inside, numpy array of booleans + """ + + tensors = self.getTensor(locType) + if type(pts) == list: + pts = np.array(pts) + assert type(pts) == np.ndarray, "must be a numpy array" + if self.dim > 1: + assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + elif len(pts.shape) == 1: + pts = pts[:,np.newaxis] + else: + assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + + inside = np.ones(pts.shape[0],dtype=bool) + for i, tensor in enumerate(tensors): + inside = inside & (pts[:,i] >= tensor.min()) & (pts[:,i] <= tensor.max()) + return inside + + def getInterpolationMat(self, loc, locType, zerosOutside=False): + """ Produces interpolation 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 interpolation matrix + + locType can be:: + + 'Ex' -> x-component of field defined on edges + 'Ey' -> y-component of field defined on edges + 'Ez' -> z-component of field defined on edges + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces + 'N' -> scalar field defined on nodes + 'CC' -> scalar field defined on cell centers + """ + + if type(loc) == list: + loc = np.array(loc) + assert type(loc) == np.ndarray, "must be a numpy array" + if self.dim > 1: + assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + elif len(loc.shape) == 1: + loc = loc[:,np.newaxis] + else: + assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + + if zerosOutside is False: + assert np.all(self.isInside(loc)), "Points outside of mesh" + else: + indZeros = np.logical_not(self.isInside(loc)) + loc[indZeros, :] = np.array([v.mean() for v in self.getTensor('CC')]) + + if locType in ['Fx','Fy','Fz','Ex','Ey','Ez']: + ind = {'x':0, 'y':1, 'z':2}[locType[1]] + assert self.dim >= ind, 'mesh is not high enough dimension.' + nF_nE = self.vnF if 'F' in locType else self.vnE + components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE] + components[ind] = Utils.interpmat(loc, *self.getTensor(locType)) + Q = sp.hstack(components) + elif locType in ['CC', 'N']: + Q = Utils.interpmat(loc, *self.getTensor(locType)) + else: + raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) + + if zerosOutside: + Q[indZeros, :] = 0 + + return Q.tocsr() + + class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): """ @@ -320,84 +401,6 @@ class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2), Utils.mkvc(l3)] return self._edge - # --------------- Methods --------------------- - - def isInside(self, pts, locType='N'): - """ - Determines if a set of points are inside a mesh. - - :param numpy.ndarray pts: Location of points to test - :rtype numpy.ndarray - :return inside, numpy array of booleans - """ - - tensors = self.getTensor(locType) - if type(pts) == list: - pts = np.array(pts) - assert type(pts) == np.ndarray, "must be a numpy array" - if self.dim > 1: - assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" - elif len(pts.shape) == 1: - pts = pts[:,np.newaxis] - else: - assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" - - inside = np.ones(pts.shape[0],dtype=bool) - for i, tensor in enumerate(tensors): - inside = inside & (pts[:,i] >= tensor.min()) & (pts[:,i] <= tensor.max()) - return inside - - def getInterpolationMat(self, loc, locType, zerosOutside=False): - """ Produces interpolation 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 interpolation matrix - - locType can be:: - - 'Ex' -> x-component of field defined on edges - 'Ey' -> y-component of field defined on edges - 'Ez' -> z-component of field defined on edges - 'Fx' -> x-component of field defined on faces - 'Fy' -> y-component of field defined on faces - 'Fz' -> z-component of field defined on faces - 'N' -> scalar field defined on nodes - 'CC' -> scalar field defined on cell centers - """ - - if type(loc) == list: - loc = np.array(loc) - assert type(loc) == np.ndarray, "must be a numpy array" - if self.dim > 1: - assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" - elif len(loc.shape) == 1: - loc = loc[:,np.newaxis] - else: - assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" - - if zerosOutside is False: - assert np.all(self.isInside(loc)), "Points outside of mesh" - else: - indZeros = np.logical_not(self.isInside(loc)) - loc[indZeros, :] = np.array([v.mean() for v in self.getTensor('CC')]) - - ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1 - if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind: - nF_nE = self.vnF if 'F' in locType else self.vnE - components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE] - components[ind] = Utils.interpmat(loc, *self.getTensor(locType)) - Q = sp.hstack(components) - elif locType in ['CC', 'N']: - Q = Utils.interpmat(loc, *self.getTensor(locType)) - else: - raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) - if zerosOutside: - Q[indZeros, :] = 0 - return Q.tocsr() - - @property def faceBoundaryInd(self): """ diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 89311bac..ed316775 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -163,7 +163,6 @@ class TestFaceDiv2D(OrderTest): F = self.M.projectFaceVector(Fc) divF = self.M.faceDiv.dot(F) - tm = Mesh.TensorMesh([self.M.nCx, self.M.nCz]) divF_anal = call3(sol, self.M.gridCC) err = np.linalg.norm((divF-divF_anal), np.inf) From 002d3a2eccc6a0aa0b1c2ba2af22ac8269f07972 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 17 Mar 2014 12:10:27 -0700 Subject: [PATCH 19/23] Tested Edge Curl --- SimPEG/Tests/test_cylMesh.py | 37 +++++++++++++++++++++++++++++------- 1 file changed, 30 insertions(+), 7 deletions(-) diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index ed316775..f1c6e423 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -141,9 +141,6 @@ call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) cyl_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] cyl_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)] cylF2 = lambda M, fx, fy: np.vstack((cyl_row2(M.gridFx, fx, fy), cyl_row2(M.gridFz, fx, fy))) -# cylE2 = lambda M, ex, ey: np.vstack((cyl_row2(M.gridEx, ex, ey), cyl_row2(M.gridEy, ex, ey))) -# cylF3 = lambda M, fx, fy, fz: np.vstack((cyl_row3(M.gridFx, fx, fy, fz), cyl_row3(M.gridFy, fx, fy, fz), cyl_row3(M.gridFz, fx, fy, fz))) -# cylE3 = lambda M, ex, ey, ez: np.vstack((cyl_row3(M.gridEx, ex, ey, ez), cyl_row3(M.gridEy, ex, ey, ez), cyl_row3(M.gridEz, ex, ey, ez))) class TestFaceDiv2D(OrderTest): @@ -177,11 +174,37 @@ class TestEdgeCurl2D(OrderTest): meshDimension = 2 def getError(self): - pass - #TODO! + # To Recreate or change the functions: - # def test_order(self): - # self.orderTest() + # import sympy + # r,t,z = sympy.symbols('r,t,z') + + # fR = 0 + # fZ = 0 + # fT = sympy.sin(2.*sympy.pi*z) + + # print 1/r*sympy.diff(fZ,t) - sympy.diff(fT,z) + # print sympy.diff(fR,z) - sympy.diff(fZ,r) + # print 1/r*(sympy.diff(r*fT,r) - sympy.diff(fR,t)) + + funT = lambda r, t, z: np.sin(2.*np.pi*z) + + solR = lambda r, z: -2.0*np.pi*np.cos(2.0*np.pi*z) + solZ = lambda r, z: np.sin(2.0*np.pi*z)/r + + E = call3(funT, self.M.gridEy) + + curlE = self.M.edgeCurl.dot(E) + + Fc = cylF2(self.M, solR, solZ) + Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]] + curlE_anal = self.M.projectFaceVector(Fc) + + err = np.linalg.norm((curlE-curlE_anal), np.inf) + return err + + def test_order(self): + self.orderTest() class TestCyl3DMesh(unittest.TestCase): From a5faa8b8919ed1b0d5017680e338e1e2707e0085 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 1 Apr 2014 09:43:14 -0700 Subject: [PATCH 20/23] Added nD to Survey.BaseRx --- SimPEG/Survey.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 936560d6..9ed4a157 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -201,6 +201,11 @@ class BaseRx(object): assert value in known, "rxType must be in ['%s']" % ("', '".join(known)) self._rxType = value + @property + def nD(self): + return self.locs.shape[0] + + class BaseTx(object): """SimPEG Transmitter Object""" From 37e957ad07478a4aa145fc8c14a7c6bb146d918d Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 14 Apr 2014 09:46:05 -0700 Subject: [PATCH 21/23] Resolve conflict.. (whoops..) --- SimPEG/Survey.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 09a2ca6c..214fbf85 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -209,10 +209,7 @@ class BaseRx(object): def nD(self): return self.locs.shape[0] -<<<<<<< HEAD -======= ->>>>>>> dbd1334e0bf48dedc12f744841e71725a9d98d50 class BaseTx(object): """SimPEG Transmitter Object""" From da7bab60b281fa0b7c2090e2e64084c33f131c73 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 14 Apr 2014 21:28:52 -0700 Subject: [PATCH 22/23] Updates to CylMesh. Innerproducts. Remove getMass. --- SimPEG/Mesh/CylMesh.py | 260 +++++++---------------------------- SimPEG/Mesh/DiffOperators.py | 60 +------- SimPEG/Mesh/TensorMesh.py | 216 +++++++++++++++-------------- SimPEG/Solver.py | 2 +- SimPEG/Tests/test_Solver.py | 2 +- SimPEG/Tests/test_cylMesh.py | 56 ++++++++ 6 files changed, 226 insertions(+), 370 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index c44afebe..f797c1c5 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -3,39 +3,10 @@ import scipy.sparse as sp from scipy.constants import pi from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, ddx, av, avExtrap from TensorMesh import BaseTensorMesh +from InnerProducts import InnerProducts - # 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 - -class CylMesh(BaseTensorMesh): +class CylMesh(BaseTensorMesh, InnerProducts): """ CylMesh is a mesh class for cylindrical problems """ @@ -49,6 +20,10 @@ class CylMesh(BaseTensorMesh): 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" + @property + def isSymmetric(self): + return self.nCy == 1 + @property def nNx(self): """ @@ -57,7 +32,7 @@ class CylMesh(BaseTensorMesh): :rtype: int :return: nNx """ - if self.nCy == 1: return self.nCx + if self.isSymmetric: return self.nCx return self.nCx + 1 @property @@ -68,7 +43,7 @@ class CylMesh(BaseTensorMesh): :rtype: int :return: nNy """ - if self.nCy == 1: return 0 + if self.isSymmetric: return 0 return self.nCy @property @@ -89,7 +64,7 @@ class CylMesh(BaseTensorMesh): :rtype: numpy.array (dim, ) :return: vnEy or None if dim < 2 """ - nNx = self.nNx if self.nCy == 1 else self.nNx - 1 + nNx = self.nNx if self.isSymmetric else self.nNx - 1 return np.r_[nNx, self.nCy, self.nNz] @property @@ -100,7 +75,7 @@ class CylMesh(BaseTensorMesh): :rtype: numpy.array (dim, ) :return: vnEz or None if nCy > 1 """ - if self.nCy == 1: + if self.isSymmetric: return np.r_[self.nNx, self.nNy, self.nCz] else: return None @@ -113,7 +88,7 @@ class CylMesh(BaseTensorMesh): :rtype: int :return: nEz """ - if self.nCy == 1: + if self.isSymmetric: return self.vnEz.prod() return (np.r_[self.nNx-1, self.nNy, self.nCz]).prod() + self.nCz @@ -130,14 +105,14 @@ class CylMesh(BaseTensorMesh): @property def vectorNx(self): """Nodal grid vector (1D) in the x direction.""" - if self.nCy == 1: + if self.isSymmetric: return self.hx.cumsum() return np.r_[0, self.hx].cumsum() @property def vectorNy(self): """Nodal grid vector (1D) in the y direction.""" - if self.nCy == 1: + if self.isSymmetric: # There aren't really any nodes, but all the grids need # somewhere to live, why not zero?! return np.r_[0] @@ -147,7 +122,7 @@ class CylMesh(BaseTensorMesh): def edge(self): """Edge lengths""" if getattr(self, '_edge', None) is None: - if self.nCy == 1: + if self.isSymmetric: self._edge = 2*pi*self.gridN[:,0] else: raise NotImplementedError('edges not yet implemented for 3D cyl mesh') @@ -186,7 +161,7 @@ class CylMesh(BaseTensorMesh): # Compute faceDivergence operator on faces D1 = self.faceDivx D3 = self.faceDivz - if self.nCy == 1: + if self.isSymmetric: D = sp.hstack((D1, D3), format="csr") elif self.nCy > 1: D2 = self.faceDivy @@ -237,7 +212,7 @@ class CylMesh(BaseTensorMesh): def nodalGrad(self): """Construct gradient operator (nodes to edges).""" # Nodal grad does not make sense for cylindrically symmetric mesh. - if self.nCy == 1: return None + if self.isSymmetric: return None raise NotImplementedError('nodalGrad not yet implemented') @property @@ -263,178 +238,49 @@ class CylMesh(BaseTensorMesh): 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 centers""" + # if getattr(self, '_aveE2CC', None) is None: + # if self.isSymmetric: + # 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 = (0.5)*sp.kron(az, ar).T + # else: + # raise NotImplementedError('wrapping in the averaging is not yet implemented') + # return self._aveE2CC + @property def aveE2CC(self): - """Averaging operator from cell edges to cell centres""" + "Construct the averaging operator on cell edges to cell centers." 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 + # The number of cell centers in each direction + n = self.vnC + if self.isSymmetric: + avR = av(n[0])[:,1:] + avR[0,0] = 1. + self._aveE2CC = (0.5)*sp.kron(av(n[2]), avR, format="csr") + else: + raise NotImplementedError('wrapping in the averaging is not yet implemented') + # self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), + # kron3(av(n[2]), speye(n[1]), av(n[0])), + # kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") return self._aveE2CC + @property def aveF2CC(self): - """Averaging operator from cell faces to cell centres""" + "Construct the averaging operator on cell faces to cell centers." 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 + n = self.vnC + if self.isSymmetric: + avR = av(n[0])[:,1:] + avR[0,0] = 1. + self._aveF2CC = (0.5)*sp.hstack((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') + # self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])), + # kron3(speye(n[2]), av(n[1]), speye(n[0])), + # kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr") 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.vectorNx.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() - - diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 15d8ed2b..f8d448d4 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -623,11 +623,11 @@ class DiffOperators(object): raise Exception('Edge Averaging does not make sense in 1D: Use Identity?') elif(self.dim == 2): self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])), - sp.kron(speye(n[1]), av(n[0]))), format="csr") + sp.kron(speye(n[1]), av(n[0]))), format="csr") elif(self.dim == 3): self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])), - kron3(av(n[2]), speye(n[1]), av(n[0])), - kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") + kron3(av(n[2]), speye(n[1]), av(n[0])), + kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr") return self._aveE2CC @property @@ -696,57 +696,3 @@ class DiffOperators(object): kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr") return self._aveN2F - # --------------- Methods --------------------- - - def getMass(self, materialProp=None, loc='e'): - """ Produces mass matricies. - - :param str loc: Average to location: 'e'-edges, 'f'-faces - :param None,float,numpy.ndarray materialProp: property to be averaged (see below) - :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*self.nCy) - 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 getFaceMassDeriv(self): - Av = self.aveF2CC - return Av.T * sdiag(self.vol) - - def getEdgeMassDeriv(self): - Av = self.aveE2CC - return Av.T * sdiag(self.vol) diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index cc3bcbb1..4d1b16a2 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -239,6 +239,8 @@ class BaseTensorMesh(BaseRectangularMesh): nF_nE = self.vnF if 'F' in locType else self.vnE components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE] components[ind] = Utils.interpmat(loc, *self.getTensor(locType)) + # remove any zero blocks (hstack complains) + components = [comp for comp in components if comp.shape[1] > 0] Q = sp.hstack(components) elif locType in ['CC', 'N']: Q = Utils.interpmat(loc, *self.getTensor(locType)) @@ -251,6 +253,116 @@ class BaseTensorMesh(BaseRectangularMesh): return Q.tocsr() + def _fastFaceInnerProduct(self, materialProperty=None, invertProperty=False): + """ + Fast version of getFaceInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :param bool invertProperty: inverts the material property + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + return self._fastInnerProduct('F', materialProperty=materialProperty, invertProperty=invertProperty) + + + def _fastEdgeInnerProduct(self, materialProperty=None, invertProperty=False): + """ + Fast version of getEdgeInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :param bool invertProperty: inverts the material property + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nE, nE) + """ + return self._fastInnerProduct('E', materialProperty=materialProperty, invertProperty=invertProperty) + + + def _fastInnerProduct(self, AvType, materialProperty=None, invertProperty=False): + """ + Fast version of getFaceInnerProduct. + This does not handle the case of a full tensor materialProperty. + + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param str AvType: 'E' or 'F' + :param bool returnP: returns the projection matrices + :param bool invertProperty: inverts the material property + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + if materialProperty is None: + materialProperty = np.ones(self.nC) + + if invertProperty: + materialProperty = 1./materialProperty + + if Utils.isScalar(materialProperty): + materialProperty = materialProperty*np.ones(self.nC) + + if materialProperty.size == self.nC: + Av = getattr(self, 'ave'+AvType+'2CC') + Vprop = self.vol * Utils.mkvc(materialProperty) + return self.dim * Utils.sdiag(Av.T * Vprop) + if materialProperty.size == self.nC*self.dim: + Av = getattr(self, 'ave'+AvType+'2CCV') + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty)) + + + def _fastFaceInnerProductDeriv(self, materialProperty=None, v=None): + """ + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + return self._fastInnerProductDeriv('F', materialProperty=materialProperty, v=v) + + + def _fastEdgeInnerProductDeriv(self, materialProperty=None, v=None): + """ + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nE, nE) + """ + return self._fastInnerProductDeriv('E', materialProperty=materialProperty, v=v) + + + def _fastInnerProductDeriv(self, AvType, materialProperty=None, v=None): + """ + :param str AvType: 'E' or 'F' + :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + if materialProperty is None: + return 0.0 + + if Utils.isScalar(materialProperty): + Av = getattr(self, 'ave'+AvType+'2CC') + V = Utils.sdiag(self.vol) + ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1)) + if v is None: + return self.dim * Av.T * V * ones + return Utils.sdiag(v) * self.dim * Av.T * V * ones + + if materialProperty.size == self.nC: + Av = getattr(self, 'ave'+AvType+'2CC') + V = Utils.sdiag(self.vol) + if v is None: + return self.dim * Av.T * V + return Utils.sdiag(v) * self.dim * Av.T * V + + if materialProperty.size == self.nC*self.dim: # anisotropic + Av = getattr(self, 'ave'+AvType+'2CCV') + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + if v is None: + return Av.T * V + return Utils.sdiag(v) * Av.T * V + + class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): """ @@ -449,111 +561,7 @@ class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): indzu = (self.gridCC[:,2]==max(self.gridCC[:,2])) return indxd, indxu, indyd, indyu, indzd, indzu - def _fastFaceInnerProduct(self, materialProperty=None, invertProperty=False): - """ - Fast version of getFaceInnerProduct. - This does not handle the case of a full tensor materialProperty. - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param bool returnP: returns the projection matrices - :param bool invertProperty: inverts the material property - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ - return self._fastInnerProduct('F', materialProperty=materialProperty, invertProperty=invertProperty) - - - def _fastEdgeInnerProduct(self, materialProperty=None, invertProperty=False): - """ - Fast version of getEdgeInnerProduct. - This does not handle the case of a full tensor materialProperty. - - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param bool returnP: returns the projection matrices - :param bool invertProperty: inverts the material property - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nE, nE) - """ - return self._fastInnerProduct('E', materialProperty=materialProperty, invertProperty=invertProperty) - - - def _fastInnerProduct(self, AvType, materialProperty=None, invertProperty=False): - """ - Fast version of getFaceInnerProduct. - This does not handle the case of a full tensor materialProperty. - - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param str AvType: 'E' or 'F' - :param bool returnP: returns the projection matrices - :param bool invertProperty: inverts the material property - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ - if materialProperty is None: - materialProperty = np.ones(self.nC) - - if invertProperty: - materialProperty = 1./materialProperty - - if Utils.isScalar(materialProperty): - materialProperty = materialProperty*np.ones(self.nC) - - if materialProperty.size == self.nC: - Av = getattr(self, 'ave'+AvType+'2CC') - Vprop = self.vol * Utils.mkvc(materialProperty) - return self.dim * Utils.sdiag(Av.T * Vprop) - if materialProperty.size == self.nC*self.dim: - Av = getattr(self, 'ave'+AvType+'2CCV') - V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty)) - - - def _fastFaceInnerProductDeriv(self, materialProperty=None, v=None): - """ - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ - return self._fastInnerProductDeriv('F', materialProperty=materialProperty, v=v) - - - def _fastEdgeInnerProductDeriv(self, materialProperty=None, v=None): - """ - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nE, nE) - """ - return self._fastInnerProductDeriv('E', materialProperty=materialProperty, v=v) - - - def _fastInnerProductDeriv(self, AvType, materialProperty=None, v=None): - """ - :param str AvType: 'E' or 'F' - :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ - if materialProperty is None: - return None - if Utils.isScalar(materialProperty): - Av = getattr(self, 'ave'+AvType+'2CC') - V = Utils.sdiag(self.vol) - ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1)) - if v is None: - return self.dim * Av.T * V * ones - return Utils.sdiag(v) * self.dim * Av.T * V * ones - if materialProperty.size == self.nC: - Av = getattr(self, 'ave'+AvType+'2CC') - V = Utils.sdiag(self.vol) - if v is None: - return self.dim * Av.T * V - return Utils.sdiag(v) * self.dim * Av.T * V - if materialProperty.size == self.nC*self.dim: # anisotropic - Av = getattr(self, 'ave'+AvType+'2CCV') - V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - if v is None: - return Av.T * V - return Utils.sdiag(v) * Av.T * V if __name__ == '__main__': diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index cc6d6338..85a87bdc 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -337,7 +337,7 @@ if __name__ == '__main__': D = M.faceDiv G = M.cellGrad - Msig = M.getFaceMass() + Msig = M.getFaceInnerProduct() A = D*Msig*G A[0,0] *= 10 # remove the constant null space from the matrix diff --git a/SimPEG/Tests/test_Solver.py b/SimPEG/Tests/test_Solver.py index f8a8983f..76177990 100644 --- a/SimPEG/Tests/test_Solver.py +++ b/SimPEG/Tests/test_Solver.py @@ -22,7 +22,7 @@ class TestSolver(unittest.TestCase): D = M.faceDiv G = M.cellGrad - Msig = M.getFaceMass() + Msig = M.getFaceInnerProduct() A = D*Msig*G A[0,0] *= 10 # remove the constant null space from the matrix diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index f1c6e423..5eaef428 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -206,6 +206,62 @@ class TestEdgeCurl2D(OrderTest): def test_order(self): self.orderTest() + +# class TestInnerProducts2D(OrderTest): +# """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" + +# meshTypes = MESHTYPES +# meshDimension = 2 +# meshSizes = [4, 8, 16, 32, 64, 128] + +# def getError(self): + +# funR = lambda r, t, z: np.cos(2.0*np.pi*z) +# funT = lambda r, t, z: 0*t +# funZ = lambda r, t, z: np.sin(2.0*np.pi*r) + +# call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) + +# sigma1 = lambda r, t, z: z+1 +# sigma2 = lambda r, t, z: r*z+50 +# sigma3 = lambda r, t, z: 3+t*r +# sigma4 = lambda r, t, z: 0.1*r*t*z +# sigma5 = lambda r, t, z: 0.2*z*r*t +# sigma6 = lambda r, t, z: 0.1*t + +# Gc = self.M.gridCC +# if self.sigmaTest == 1: +# sigma = np.c_[call(sigma1, Gc)] +# analytic = 144877./360 # Found using sympy. z=5 +# elif self.sigmaTest == 2: +# sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)] +# analytic = 189959./120 # Found using sympy. z=5 +# elif self.sigmaTest == 3: +# sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)] +# analytic = 781427./360 # Found using sympy. z=5 + +# if self.location == 'edges': +# E = call(funT, self.M.gridEy) +# A = self.M.getEdgeInnerProduct(sigma) +# numeric = E.T.dot(A.dot(E)) +# elif self.location == 'faces': +# Fr = call(funR, self.M.gridFx) +# Fz = call(funZ, self.M.gridFz) +# A = self.M.getFaceInnerProduct(sigma) +# F = np.r_[Fr,Fz] +# numeric = F.T.dot(A.dot(F)) + +# print numeric +# err = np.abs(numeric - analytic) +# return err + +# def test_order1_faces(self): +# self.name = "2D Face Inner Product - Isotropic" +# self.location = 'faces' +# self.sigmaTest = 1 +# self.orderTest() + + class TestCyl3DMesh(unittest.TestCase): def setUp(self): From bb1be7b96ecb549a958bdcac4eac03c0ef98a24a Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 15 Apr 2014 09:05:07 -0700 Subject: [PATCH 23/23] Added Cyl1DMesh back in for depreciation purposes. --- SimPEG/Mesh/Cyl1DMesh.py | 477 +++++++++++++++++++++++++++++++++++++++ SimPEG/Mesh/CylMesh.py | 7 + SimPEG/Mesh/__init__.py | 1 + 3 files changed, 485 insertions(+) create mode 100644 SimPEG/Mesh/Cyl1DMesh.py diff --git a/SimPEG/Mesh/Cyl1DMesh.py b/SimPEG/Mesh/Cyl1DMesh.py new file mode 100644 index 00000000..d53a06bf --- /dev/null +++ b/SimPEG/Mesh/Cyl1DMesh.py @@ -0,0 +1,477 @@ +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" + print 'PendingDeprecationWarning: Cyl1DMesh has been replaced by CylMesh. hy == theta == 2pi, use one cell to make this cylindrically symmetric.' + 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 index f797c1c5..f9fb0503 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -9,6 +9,13 @@ from InnerProducts import InnerProducts class CylMesh(BaseTensorMesh, InnerProducts): """ CylMesh is a mesh class for cylindrical problems + + :: + + cs, nc, npad = 20., 30, 8 + hx = Utils.meshTensors(((npad+10,cs,0.7), (nc,cs), (npad,cs))) + hz = Utils.meshTensors(((npad,cs), (nc,cs), (npad,cs))) + mesh = Mesh.CylMesh([hx,1,hz], [0.,0,-hz.sum()/2.]) """ _meshType = 'CYL' diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py index 38eaf527..a79e93dc 100644 --- a/SimPEG/Mesh/__init__.py +++ b/SimPEG/Mesh/__init__.py @@ -1,4 +1,5 @@ from TensorMesh import TensorMesh from CylMesh import CylMesh +from Cyl1DMesh import Cyl1DMesh from LogicallyRectMesh import LogicallyRectMesh from TreeMesh import TreeMesh