diff --git a/AUTHORS.rst b/AUTHORS.rst new file mode 100644 index 00000000..221c85c2 --- /dev/null +++ b/AUTHORS.rst @@ -0,0 +1,8 @@ +Credits +======= + +- Rowan Cockett +- Eldad Haber +- Lindsey Heagy +- Seogi Kang +- Dave Marchant diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py index a1b631d1..9083e59d 100644 --- a/SimPEG/Mesh/BaseMesh.py +++ b/SimPEG/Mesh/BaseMesh.py @@ -26,7 +26,7 @@ class BaseMesh(object): # Ensure x0 & n are 1D vectors self._n = np.array(n, dtype=int).ravel() - self._x0 = np.array(x0).ravel() + self._x0 = np.array(x0, dtype=float).ravel() @property def x0(self): @@ -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 @@ -403,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/Cyl1DMesh.py b/SimPEG/Mesh/Cyl1DMesh.py index f27d3d0f..d53a06bf 100644 --- a/SimPEG/Mesh/Cyl1DMesh.py +++ b/SimPEG/Mesh/Cyl1DMesh.py @@ -12,6 +12,7 @@ class Cyl1DMesh(object): 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" diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py new file mode 100644 index 00000000..f9fb0503 --- /dev/null +++ b/SimPEG/Mesh/CylMesh.py @@ -0,0 +1,293 @@ +import numpy as np +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 + + +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' + + _unitDimensions = [1, 2*np.pi, 1] + + def __init__(self, h, x0=None): + BaseTensorMesh.__init__(self, h, x0) + assert self.dim == 3, "dim of mesh must equal 3, for a cylindrically symmetric mesh use [hx, 1, hz]" + assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi" + + @property + def isSymmetric(self): + return self.nCy == 1 + + @property + def nNx(self): + """ + Number of nodes in the x-direction + + :rtype: int + :return: nNx + """ + if self.isSymmetric: return self.nCx + return self.nCx + 1 + + @property + def nNy(self): + """ + Number of nodes in the y-direction + + :rtype: int + :return: nNy + """ + if self.isSymmetric: return 0 + return self.nCy + + @property + def vnFx(self): + """ + Number of x-faces in each direction + + :rtype: numpy.array (dim, ) + :return: vnFx + """ + return self.vnC + + @property + def vnEy(self): + """ + Number of y-edges in each direction + + :rtype: numpy.array (dim, ) + :return: vnEy or None if dim < 2 + """ + nNx = self.nNx if self.isSymmetric 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.isSymmetric: + 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.isSymmetric: + 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.""" + 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 x direction.""" + 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.isSymmetric: + # 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 + + @property + def edge(self): + """Edge lengths""" + if getattr(self, '_edge', None) is None: + if self.isSymmetric: + self._edge = 2*pi*self.gridN[:,0] + else: + 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: + 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 + + @property + def vol(self): + """Volume of each cell""" + if getattr(self, '_vol', None) is None: + 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 + + #################################################### + # Operators + #################################################### + + @property + def faceDiv(self): + """Construct divergence operator (face-stg to cell-centres).""" + if getattr(self, '_faceDiv', None) is None: + n = self.vnC + # Compute faceDivergence operator on faces + D1 = self.faceDivx + D3 = self.faceDivz + if self.isSymmetric: + 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).""" + 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('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).""" + 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 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).""" + # Nodal grad does not make sense for cylindrically symmetric mesh. + if self.isSymmetric: return None + 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.""" + 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") + 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 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): + "Construct the averaging operator on cell edges to cell centers." + if getattr(self, '_aveE2CC', None) is None: + # 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): + "Construct the averaging operator on cell faces to cell centers." + if getattr(self, '_aveF2CC', None) is None: + 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 diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index a7a7b0cb..f8d448d4 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) @@ -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/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 22c6a521..44d5f37d 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -284,7 +284,6 @@ class InnerProducts(object): # | |/ # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) - def _getFacePx(M): assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh' return _getFacePx_Rectangular(M) diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 1bc7cadb..4d1b16a2 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -1,50 +1,31 @@ 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 -class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): - """ - TensorMesh is a mesh class that deals with tensor product meshes. +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) - 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' assert len(h_in) in [1,2,3], 'h_in must be of dimension 1, 2, or 3' 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 = 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. @@ -55,279 +36,101 @@ 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 + @property + def h(self): + """h is a list containing the cell widths of the tensor mesh in each dimension.""" + return self._h - if n == 1: - outStr = outStr + ' {0:.2f},'.format(h) - else: - outStr = outStr + ' {0:d}*{1:.2f},'.format(n,h) + @property + def hx(self): + "Width of cells in the x direction" + return self._h[0] - return outStr[:-1] + @property + def hy(self): + "Width of cells in the y direction" + return None if self.dim < 2 else self._h[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:') + @property + def hz(self): + "Width of cells in the z direction" + return None if self.dim < 3 else self._h[2] - return outStr + @property + def vectorNx(self): + """Nodal grid vector (1D) in the x direction.""" + return np.r_[0., self.hx.cumsum()] + self.x0[0] - 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 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 hx(): - doc = "Width of cells in the x direction" - fget = lambda self: self._h[0] - return locals() - hx = property(**hx()) + @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 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 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 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 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 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 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 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 gridCC(self): + """Cell-centered grid.""" + return _getTensorGrid(self, 'CC') - 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 gridN(self): + """Nodal grid.""" + return _getTensorGrid(self, 'N') - 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 gridFx(self): + """Face staggered grid in the x direction.""" + if self.nFx == 0: return + return _getTensorGrid(self, 'Fx') - 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 gridFy(self): + """Face staggered grid in the y direction.""" + if self.nFy == 0 or self.dim < 2: return + return _getTensorGrid(self, 'Fy') - 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 gridFz(self): + """Face staggered grid in the z direction.""" + if self.nFz == 0 or self.dim < 3: return + return _getTensorGrid(self, 'Fz') - def gridCC(): - doc = "Cell-centered grid." + @property + def gridEx(self): + """Edge staggered grid in the x direction.""" + if self.nEx == 0: return + return _getTensorGrid(self, 'Ex') - 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 gridEy(self): + """Edge staggered grid in the y direction.""" + if self.nEy == 0 or self.dim < 2: return + return _getTensorGrid(self, 'Ey') - def gridN(): - doc = "Nodal grid." - - 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()) - - def gridFx(): - doc = "Face staggered grid in the x direction." - - 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()) - - def gridFy(): - doc = "Face staggered grid in the y direction." - - 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()) - - # --------------- Geometries --------------------- - def vol(): - doc = "Construct cell volumes of the 3D model as 1d array." - - 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()) - - 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()) - - # --------------- Methods --------------------- + @property + def gridEz(self): + """Edge staggered grid in the z direction.""" + if self.nEz == 0 or self.dim < 3: return + return _getTensorGrid(self, 'Ez') def getTensor(self, locType): """ Returns a tensor list. @@ -367,6 +170,7 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): return [t for t in ten if t is not None] + # --------------- Methods --------------------- def isInside(self, pts, locType='N'): """ @@ -429,69 +233,26 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): 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: + 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)) + # 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)) else: raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) + if zerosOutside: Q[indZeros, :] = 0 + return Q.tocsr() - @property - def faceBoundaryInd(self): - """ - Find indices of boundary faces in each direction - """ - if self.dim==1: - indxd = (self.gridFx==min(self.gridFx)) - indxu = (self.gridFx==max(self.gridFx)) - return indxd, indxu - elif self.dim==2: - indxd = (self.gridFx[:,0]==min(self.gridFx[:,0])) - indxu = (self.gridFx[:,0]==max(self.gridFx[:,0])) - indyd = (self.gridFy[:,1]==min(self.gridFy[:,1])) - indyu = (self.gridFy[:,1]==max(self.gridFy[:,1])) - return indxd, indxu, indyd, indyu - elif self.dim==3: - indxd = (self.gridFx[:,0]==min(self.gridFx[:,0])) - indxu = (self.gridFx[:,0]==max(self.gridFx[:,0])) - indyd = (self.gridFy[:,1]==min(self.gridFy[:,1])) - indyu = (self.gridFy[:,1]==max(self.gridFy[:,1])) - indzd = (self.gridFz[:,2]==min(self.gridFz[:,2])) - indzu = (self.gridFz[:,2]==max(self.gridFz[:,2])) - return indxd, indxu, indyd, indyu, indzd, indzu - - @property - def cellBoundaryInd(self): - """ - Find indices of boundary faces in each direction - """ - if self.dim==1: - indxd = (self.gridCC==min(self.gridCC)) - indxu = (self.gridCC==max(self.gridCC)) - return indxd, indxu - elif self.dim==2: - indxd = (self.gridCC[:,0]==min(self.gridCC[:,0])) - indxu = (self.gridCC[:,0]==max(self.gridCC[:,0])) - indyd = (self.gridCC[:,1]==min(self.gridCC[:,1])) - indyu = (self.gridCC[:,1]==max(self.gridCC[:,1])) - return indxd, indxu, indyd, indyu - elif self.dim==3: - indxd = (self.gridCC[:,0]==min(self.gridCC[:,0])) - indxu = (self.gridCC[:,0]==max(self.gridCC[:,0])) - indyd = (self.gridCC[:,1]==min(self.gridCC[:,1])) - indyu = (self.gridCC[:,1]==max(self.gridCC[:,1])) - indzd = (self.gridCC[:,2]==min(self.gridCC[:,2])) - 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. @@ -577,7 +338,8 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): :return: M, the inner product matrix (nF, nF) """ if materialProperty is None: - return None + return 0.0 + if Utils.isScalar(materialProperty): Av = getattr(self, 'ave'+AvType+'2CC') V = Utils.sdiag(self.vol) @@ -585,12 +347,14 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): 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)) @@ -599,6 +363,207 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): return Utils.sdiag(v) * Av.T * V + +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): + """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 + + @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 + + @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 + + @property + def faceBoundaryInd(self): + """ + Find indices of boundary faces in each direction + """ + if self.dim==1: + indxd = (self.gridFx==min(self.gridFx)) + indxu = (self.gridFx==max(self.gridFx)) + return indxd, indxu + elif self.dim==2: + indxd = (self.gridFx[:,0]==min(self.gridFx[:,0])) + indxu = (self.gridFx[:,0]==max(self.gridFx[:,0])) + indyd = (self.gridFy[:,1]==min(self.gridFy[:,1])) + indyu = (self.gridFy[:,1]==max(self.gridFy[:,1])) + return indxd, indxu, indyd, indyu + elif self.dim==3: + indxd = (self.gridFx[:,0]==min(self.gridFx[:,0])) + indxu = (self.gridFx[:,0]==max(self.gridFx[:,0])) + indyd = (self.gridFy[:,1]==min(self.gridFy[:,1])) + indyu = (self.gridFy[:,1]==max(self.gridFy[:,1])) + indzd = (self.gridFz[:,2]==min(self.gridFz[:,2])) + indzu = (self.gridFz[:,2]==max(self.gridFz[:,2])) + return indxd, indxu, indyd, indyu, indzd, indzu + + @property + def cellBoundaryInd(self): + """ + Find indices of boundary faces in each direction + """ + if self.dim==1: + indxd = (self.gridCC==min(self.gridCC)) + indxu = (self.gridCC==max(self.gridCC)) + return indxd, indxu + elif self.dim==2: + indxd = (self.gridCC[:,0]==min(self.gridCC[:,0])) + indxu = (self.gridCC[:,0]==max(self.gridCC[:,0])) + indyd = (self.gridCC[:,1]==min(self.gridCC[:,1])) + indyu = (self.gridCC[:,1]==max(self.gridCC[:,1])) + return indxd, indxu, indyd, indyu + elif self.dim==3: + indxd = (self.gridCC[:,0]==min(self.gridCC[:,0])) + indxu = (self.gridCC[:,0]==max(self.gridCC[:,0])) + indyd = (self.gridCC[:,1]==min(self.gridCC[:,1])) + indyu = (self.gridCC[:,1]==max(self.gridCC[:,1])) + indzd = (self.gridCC[:,2]==min(self.gridCC[:,2])) + indzu = (self.gridCC[:,2]==max(self.gridCC[:,2])) + return indxd, indxu, indyd, indyu, indzd, indzu + + + + if __name__ == '__main__': print('Welcome to tensor mesh!') diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/View.py similarity index 83% rename from SimPEG/Mesh/TensorView.py rename to SimPEG/Mesh/View.py index 025e27ec..d90475f5 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/View.py @@ -516,6 +516,106 @@ 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() + + if __name__ == '__main__': from SimPEG import * mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2))) @@ -525,6 +625,7 @@ if __name__ == '__main__': q = Utils.mkvc(q) A = M.faceDiv*M.cellGrad b = Solver(A).solve(q) + M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, pcolorOpts={'alpha':0.8}) M2 = Mesh.TensorMesh([10,20],x0=[10,5]) f = np.r_[np.sin(M2.gridFx[:,0]*2*np.pi), np.sin(M2.gridFy[:,1]*2*np.pi)] diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py index 0fa60201..a79e93dc 100644 --- a/SimPEG/Mesh/__init__.py +++ b/SimPEG/Mesh/__init__.py @@ -1,8 +1,5 @@ -from Cyl1DMesh import Cyl1DMesh from TensorMesh import TensorMesh -from TreeMesh import TreeMesh +from CylMesh import CylMesh +from Cyl1DMesh import Cyl1DMesh from LogicallyRectMesh import LogicallyRectMesh -from BaseMesh import BaseMesh, BaseRectangularMesh -from TensorView import TensorView -from InnerProducts import InnerProducts -from DiffOperators import DiffOperators +from TreeMesh import TreeMesh 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/Survey.py b/SimPEG/Survey.py index 09c45083..214fbf85 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -209,6 +209,7 @@ class BaseRx(object): def nD(self): return self.locs.shape[0] + class BaseTx(object): """SimPEG Transmitter Object""" diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index c4d11443..0fb1d200 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, LogicallyRectMesh +from SimPEG.Mesh import TensorMesh, LogicallyRectMesh, 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)*nc*0.5 + nc*0.5 h2 = np.random.rand(nc)*nc*0.5 + nc*0.5 @@ -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 'LRM' in self._meshType: if 'uniform' in self._meshType: kwrd = 'rect' 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_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 new file mode 100644 index 00000000..5eaef428 --- /dev/null +++ b/SimPEG/Tests/test_cylMesh.py @@ -0,0 +1,327 @@ +import unittest +import sys +from SimPEG import * +from TestUtils import OrderTest + + +class TestCyl2DMesh(unittest.TestCase): + + def setUp(self): + hx = np.r_[1,1,0.5] + hz = np.r_[2,1] + self.mesh = Mesh.CylMesh([hx, 1,hz]) + + 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) + 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])) + + 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) + 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) + 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) + 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_[1, 2, 2.5] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNx)) == 0) + 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_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] + 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) + + 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) + + def test_gridSizes(self): + self.assertTrue(self.mesh.gridCC.shape == (self.mesh.nC, 3)) + self.assertTrue(self.mesh.gridN.shape == (9, 3)) + + self.assertTrue(self.mesh.gridFx.shape == (self.mesh.nFx, 3)) + self.assertTrue(self.mesh.gridFy is None) + self.assertTrue(self.mesh.gridFz.shape == (self.mesh.nFz, 3)) + + self.assertTrue(self.mesh.gridEx is None) + self.assertTrue(self.mesh.gridEy.shape == (self.mesh.nEy, 3)) + self.assertTrue(self.mesh.gridEz is None) + + def test_gridCC(self): + 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] + 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) + 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.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] + 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) + + def test_lightOperators(self): + self.assertTrue(self.mesh.nodalGrad is None) + + + +MESHTYPES = ['uniformCylMesh'] +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)] +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))) + + +class TestFaceDiv2D(OrderTest): + name = "FaceDiv" + meshTypes = MESHTYPES + meshDimension = 2 + + def getError(self): + + funR = lambda r, z: np.sin(2.*np.pi*r) + funZ = lambda r, z: np.sin(2.*np.pi*z) + + 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) + + Fc = cylF2(self.M, funR, funZ) + Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]] + F = self.M.projectFaceVector(Fc) + + divF = self.M.faceDiv.dot(F) + divF_anal = call3(sol, self.M.gridCC) + + err = np.linalg.norm((divF-divF_anal), np.inf) + return err + + def test_order(self): + self.orderTest() + +class TestEdgeCurl2D(OrderTest): + name = "EdgeCurl" + meshTypes = MESHTYPES + meshDimension = 2 + + def getError(self): + # To Recreate or change the functions: + + # 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 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): + 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_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 == 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] + 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__': + unittest.main() diff --git a/SimPEG/Tests/test_operators.py b/SimPEG/Tests/test_operators.py index 1cf99402..1610f975 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', 'uniformLRM', 'rotateLRM'] call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 94d016f6..05796ef4 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -1,5 +1,5 @@ from matutils import * -from meshutils import exampleLrmGrid, meshTensors, points2nodes +from meshutils import exampleLrmGrid, meshTensors, points2nodes, writeUBCTensorMesh, writeUBCTensorModel from lrmutils import volTetra, faceInfo, indexCube from interputils import interpmat from ipythonutils import easyAnimate as animate diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 155a6d9f..fdbece8f 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -24,7 +24,6 @@ def exampleLrmGrid(nC, exType): amt[amt < 0] = 0 return [X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt] - def meshTensors(*args): """ **meshTensors** takes any number of tuples that have the form:: @@ -58,7 +57,7 @@ def points2nodes(mesh, pts): Move a list of the nearest nodes to a set of points :param simpeg.Mesh.TensorMesh mesh: The mesh - :param numpy.ndarray pts: Points to move} + :param numpy.ndarray pts: Points to move :rtype: numpy.ndarray :return: nodeInds """ @@ -76,6 +75,48 @@ def points2nodes(mesh, pts): return nodeInds +def writeUBCTensorMesh(mesh, fileName): + """ + Writes a SimPEG TensorMesh to a UBC-GIF format mesh file. + + :param simpeg.Mesh.TensorMesh mesh: The mesh + :param str fileName: File to write to + """ + assert mesh.dim == 3 + s = '' + s += '%i %i %i\n' %tuple(mesh.vnC) + origin = mesh.x0 + origin.dtype = float + origin[2] = origin[2]+mesh.hz.sum() + s += '%.2f %.2f %.2f\n' %tuple(origin) + s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx) + s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy) + s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz) + f = open(fileName, 'w') + f.write(s) + f.close() + +def writeUBCTensorModel(mesh, model, fileName): + """ + Writes a model associated with a SimPEG TensorMesh + to a UBC-GIF format model file. + + :param simpeg.Mesh.TensorMesh mesh: The mesh + :param numpy.ndarray model: The model + :param str fileName: File to write to + """ + + + + # Reshape to [z,y,x] + model3D = np.reshape(model, mesh.vnC) + # Permute to [z,x,y] + model3D = np.swapaxes(model3D, 1, 2) + # Flip z to positive down + model3D = model3D[::-1,:,:] + + np.savetxt(fileName, mkvc(model3D)) + if __name__ == '__main__': from SimPEG import Mesh