#66: Update MeshTensor to be supported directly from TensorMesh. Also added 'C' to center tensors easily in the x0 option.

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)

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.])
     """
 

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::

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)

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"

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

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()

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,