From c9111e4151b874d10fb8263f84f54a97c6d801ad Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 26 Apr 2014 12:16:49 -0700 Subject: [PATCH] #66: Update MeshTensor to be supported directly from TensorMesh. Also added 'C' to center tensors easily in the x0 option. --- SimPEG/Maps.py | 2 +- SimPEG/Mesh/CylMesh.py | 4 +-- SimPEG/Mesh/TensorMesh.py | 28 ++++++++++++--- SimPEG/Mesh/View.py | 12 ++++--- SimPEG/Tests/test_tensorMesh.py | 6 +++- SimPEG/Utils/__init__.py | 2 +- SimPEG/Utils/meshutils.py | 62 ++++++++++++++++++++++++--------- docs/api_Mesh.rst | 4 +-- 8 files changed, 87 insertions(+), 33 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 37b210d9..a32ff13b 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -251,7 +251,7 @@ class Mesh2Mesh(IdentityMap): from SimPEG import * M = Mesh.TensorMesh([100,100]) - h1 = Utils.meshTensors(((7,6,1.5),(10,6),(7,6,1.5))) + h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)]) h1 = h1/h1.sum() M2 = Mesh.TensorMesh([h1,h1]) V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index f9fb0503..db40a9ce 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -13,8 +13,8 @@ class CylMesh(BaseTensorMesh, InnerProducts): :: 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))) + hx = Utils.meshTensor([(cs,npad+10,-0.7), (cs,nc), (cs,npad,1.3)]) + hz = Utils.meshTensor([(cs,npad ,-1.3), (cs,nc), (cs,npad,1.3)]) mesh = Mesh.CylMesh([hx,1,hz], [0.,0,-hz.sum()/2.]) """ diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 72166b1f..c97aab82 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -18,20 +18,37 @@ class BaseTensorMesh(BaseRectangularMesh): _unitDimensions = [1, 1, 1] - def __init__(self, h_in, x0=None): + def __init__(self, h_in, x0_in=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, np.int_]: + if Utils.isScalar(h_i) and type(h_i) is not np.ndarray: # This gives you something over the unit cube. h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i) + elif type(h_i) is list: + h_i = Utils.meshTensor(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. + x0 = np.zeros(len(h)) + if x0_in is not None: + assert len(h) == len(x0_in), "Dimension mismatch. x0 != len(h)" + for i in range(len(h)): + x_i, h_i = x0_in[i], h[i] + if Utils.isScalar(x_i): + x0[i] = x_i + elif x_i == '0': + x0[i] = 0.0 + elif x_i == 'C': + x0[i] = -h_i.sum()*0.5 + elif x_i == 'N': + x0[i] = -h_i.sum() + else: + raise Exception("x0[%i] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative." % i) + BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0) - assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)" # Ensure h contains 1D vectors self._h = [Utils.mkvc(x.astype(float)) for x in h] @@ -376,12 +393,13 @@ class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): mesh = Mesh.TensorMesh([hx, hy, hz]) - Example of a padded tensor mesh: + Example of a padded tensor mesh using :func:`SimPEG.Utils.meshTensor`: .. plot:: + :include-source: from SimPEG import Mesh, Utils - M = Mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + M = Mesh.TensorMesh([(10,10,-1.3),(10,40),(10,10,1.3)], [(10,10,-1.3),(10,20)]) M.plotGrid() For a quick tensor mesh on a (10x12x15) unit cube:: diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index d90475f5..cc1c7279 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -178,8 +178,10 @@ class TensorView(object): .. plot:: from SimPEG import * - mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2))) - M = Mesh.TensorMesh(mT) + hx = [(5,2,-1.3),(2,4),(5,2,1.3)] + hy = [(2,2,-1.3),(2,6),(2,2,1.3)] + hz = [(2,2,-1.3),(2,6),(2,2,1.3)] + M = Mesh.TensorMesh([hx,hy,hz]) q = np.zeros(M.vnC) q[[4,4],[4,4],[2,6]]=[-1,1] q = Utils.mkvc(q) @@ -618,8 +620,10 @@ class LomView(object): 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))) - M = Mesh.TensorMesh(mT, x0=[10,20,14]) + hx = [(5,2,-1.3),(2,4),(5,2,1.3)] + hy = [(2,2,-1.3),(2,6),(2,2,1.3)] + hz = [(2,2,-1.3),(2,6),(2,2,1.3)] + M = Mesh.TensorMesh([hx,hy,hz], x0=[10,20,14]) q = np.zeros(M.vnC) q[[4,4],[4,4],[2,6]]=[-1,1] q = Utils.mkvc(q) diff --git a/SimPEG/Tests/test_tensorMesh.py b/SimPEG/Tests/test_tensorMesh.py index 45c03399..23658739 100644 --- a/SimPEG/Tests/test_tensorMesh.py +++ b/SimPEG/Tests/test_tensorMesh.py @@ -11,7 +11,7 @@ class BasicTensorMeshTests(unittest.TestCase): a = np.array([1, 1, 1]) b = np.array([1, 2]) c = np.array([1, 4]) - self.mesh2 = TensorMesh([a, b], np.array([3, 5])) + self.mesh2 = TensorMesh([a, b], [3, 5]) self.mesh3 = TensorMesh([a, b, c]) def test_vectorN_2D(self): @@ -55,6 +55,10 @@ class BasicTensorMeshTests(unittest.TestCase): t1 = np.all(self.mesh2.edge == test_edge) self.assertTrue(t1) + def test_oneCell(self): + hx = np.array([1e-5]) + M = TensorMesh([hx]) + self.assertTrue(M.nC == 1) class TestPoissonEqn(OrderTest): name = "Poisson Equation" diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index eb8975b9..85effc55 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -1,6 +1,6 @@ from matutils import * from codeutils import * -from meshutils import exampleLrmGrid, meshTensors, closestPoints, writeUBCTensorMesh, writeUBCTensorModel +from meshutils import exampleLrmGrid, meshTensor, closestPoints, 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 b64cecfc..43efcc33 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -2,6 +2,8 @@ import numpy as np from scipy import sparse as sp from matutils import mkvc, ndgrid, sub2ind, sdiag from codeutils import asArray_N_x_Dim +from codeutils import isScalar + def exampleLrmGrid(nC, exType): assert type(nC) == list, "nC must be a list containing the number of nodes" @@ -25,33 +27,57 @@ 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): +def meshTensor(value): """ - **meshTensors** takes any number of tuples that have the form:: + **meshTensor** takes a list of numbers and tuples that have the form:: - mT = ( (numPad, sizeStart [, increaseFactor]), (numCore, sizeCore), (numPad, sizeStart [, increaseFactor]) ) + mT = [ float, (cellSize, numCell), (cellSize, numCell, factor) ] - .. note:: + For example, a time domain mesh code needs + many time steps at one time:: - The increaseFactor is an optional input. + [(1e-5, 30), (1e-4, 30), 1e-3] + Means take 30 steps at 1e-5 and then 30 more at 1e-4, + and then one step of 1e-3. + + Tensor meshes can also be created by increase factors:: + + [(10.0, 5, -1.3), (10.0, 50), (10.0, 5, 1.3)] + + When there is a third number in the tuple, it + refers to the increase factor, if this number + is negative this section of the tensor is flipped right-to-left. .. 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() + from SimPEG import Mesh + tx = [(10.0,10,-1.3),(10.0,40),(10.0,10,1.3)] + ty = [(10.0,10,-1.3),(10.0,40)] + M = Mesh.TensorMesh([tx, ty]) + M.plotGrid(showIt=True) """ - def padding(num, start, factor=1.3, reverse=False): - pad = ((np.ones(num)*factor)**np.arange(num))*start - if reverse: pad = pad[::-1] - return pad - tensors = tuple() - for i, arg in enumerate(args): - tensors += (np.r_[padding(*arg[0],reverse=True),np.ones(arg[1][0])*arg[1][1],padding(*arg[2])],) + if type(value) is not list: + raise Exception('meshTensor must be a list of scalars and tuples.') - return list(tensors) if len(tensors) > 1 else tensors[0] + proposed = [] + for v in value: + if isScalar(v): + proposed += [float(v)] + elif type(v) is tuple and len(v) == 2: + proposed += [float(v[0])]*int(v[1]) + elif type(v) is tuple and len(v) == 3: + start = float(v[0]) + num = int(v[1]) + factor = float(v[2]) + pad = ((np.ones(num)*np.abs(factor))**(np.arange(num)+1))*start + if factor < 0: pad = pad[::-1] + proposed += pad.tolist() + else: + raise Exception('meshTensor must contain only scalars and len(2) or len(3) tuples.') + + return np.array(proposed) def closestPoints(mesh, pts, gridLoc='CC'): """ @@ -169,7 +195,9 @@ def writeUBCTensorModel(mesh, model, fileName): if __name__ == '__main__': from SimPEG import Mesh import matplotlib.pyplot as plt - M = Mesh.TensorMesh(meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + tx = [(10.0,10,-1.3),(10.0,40),(10.0,10,1.3)] + ty = [(10.0,10,-1.3),(10.0,40)] + M = Mesh.TensorMesh([tx, ty]) M.plotGrid() plt.gca().axis('tight') plt.show() diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst index 45172365..c947a3fe 100644 --- a/docs/api_Mesh.rst +++ b/docs/api_Mesh.rst @@ -185,8 +185,8 @@ notation:: :include-source: from SimPEG import Mesh, Utils - h1 = (5, 10, 1.5), (20, 5), (3, 10) - M = Mesh.TensorMesh(Utils.meshTensors(h1, h1)) + h1 = [(10, 5, -1.3), (5, 20), (10, 3, 1.3)] + M = Mesh.TensorMesh([h1, h1]) M.plotGrid(showIt=True) Hopefully, you now know how to create TensorMesh objects in SimPEG,