Merge branch 'em/dev' into em/rx

# Conflicts:
#	SimPEG/EM/FDEM/FDEM.py
#	SimPEG/Mesh/TensorMesh.py

tests/base/test_regularization.py

@@ -4,11 +4,17 @@ from SimPEG import *
 from scipy.sparse.linalg import dsolve
 import inspect
 
+TOL = 1e-20
 
 class RegularizationTests(unittest.TestCase):
 
     def setUp(self):
-        self.mesh2 = Mesh.TensorMesh([3, 2])
+        hx, hy, hz = np.random.rand(10), np.random.rand(9), np.random.rand(8)
+        hx, hy, hz = hx/hx.sum(), hy/hy.sum(), hz/hz.sum()
+        mesh1 = Mesh.TensorMesh([hx])
+        mesh2 = Mesh.TensorMesh([hx, hy])
+        mesh3 = Mesh.TensorMesh([hx, hy, hz])
+        self.meshlist = [mesh1,mesh2, mesh3]
 
     def test_regularization(self):
         for R in dir(Regularization):
@@ -16,18 +22,63 @@ class RegularizationTests(unittest.TestCase):
             if not inspect.isclass(r): continue
             if not issubclass(r, Regularization.BaseRegularization):
                 continue
-            # if 'Regularization' not in R: continue
-            mapping = r.mapPair(self.mesh2)
-            reg = r(self.mesh2, mapping=mapping)
-            m = np.random.rand(mapping.nP)
-            reg.mref = m[:]*np.mean(m)
 
-            print 'Check:', R
-            passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
-            self.assertTrue(passed)
-            print 'Check 2 Deriv:', R
-            passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
-            self.assertTrue(passed)
+            for i, mesh in enumerate(self.meshlist):
+
+                print 'Testing %iD'%mesh.dim
+
+                mapping = r.mapPair(mesh)
+                reg = r(mesh, mapping=mapping)
+                m = np.random.rand(mapping.nP)
+                reg.mref = np.ones_like(m)*np.mean(m)
+
+                print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                passed = reg.eval(reg.mref) < TOL
+                self.assertTrue(passed)
+
+                print 'Check:', R
+                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
+
+                print 'Check 2 Deriv:', R
+                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
+
+    def test_regularization_ActiveCells(self):
+        for R in dir(Regularization):
+            r = getattr(Regularization, R)
+            if not inspect.isclass(r): continue
+            if not issubclass(r, Regularization.BaseRegularization):
+                continue
+
+            for i, mesh in enumerate(self.meshlist):
+
+                print 'Testing Active Cells %iD'%(mesh.dim)
+
+                if mesh.dim == 1:
+                    indAct = Utils.mkvc(mesh.gridCC <= 0.8)
+                elif mesh.dim == 2:
+                    indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5)
+                elif mesh.dim == 3:
+                    indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
+
+                mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size)
+
+                reg = r(mesh, mapping=mapping, indActive=indAct)
+                m = np.random.rand(mesh.nC)[indAct]
+                reg.mref = np.ones_like(m)*np.mean(m)
+
+                print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                passed = reg.eval(reg.mref) < TOL
+                self.assertTrue(passed)
+
+                print 'Check:', R
+                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
+
+                print 'Check 2 Deriv:', R
+                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                self.assertTrue(passed)
 
 
 if __name__ == '__main__':

tests/mesh/test_MeshIO.py

@@ -0,0 +1,100 @@
+import numpy as np
+import unittest, os
+import SimPEG as simpeg
+from SimPEG.Mesh import TensorMesh, TreeMesh
+
+
+class TestTensorMeshIO(unittest.TestCase):
+
+    def setUp(self):
+        h = np.ones(16)
+        mesh = TensorMesh([h,2*h,3*h])
+        self.mesh = mesh
+
+    def test_UBCfiles(self):
+
+        mesh = self.mesh
+        # Make a vector
+        vec = np.arange(mesh.nC)
+        # Write and read
+        mesh.writeUBC('temp.msh', {'arange.txt':vec})
+        meshUBC = TensorMesh.readUBC('temp.msh')
+        vecUBC  = meshUBC.readModelUBC('arange.txt')
+
+        # The mesh
+        assert mesh.__str__() == meshUBC.__str__()
+        assert np.sum(mesh.gridCC - meshUBC.gridCC) == 0
+        assert np.sum(vec - vecUBC) == 0
+        assert np.all(np.array(mesh.h) - np.array(meshUBC.h) == 0)
+
+
+        vecUBC  = mesh.readModelUBC('arange.txt')
+        assert np.sum(vec - vecUBC) == 0
+
+        mesh.writeModelUBC('arange2.txt', vec + 1)
+        vec2UBC  = mesh.readModelUBC('arange2.txt')
+        assert np.sum(vec + 1 - vec2UBC) == 0
+
+        print 'IO of UBC tensor mesh files is working'
+        os.remove('temp.msh')
+        os.remove('arange.txt')
+        os.remove('arange2.txt')
+
+    def test_VTKfiles(self):
+        mesh = self.mesh
+        vec = np.arange(mesh.nC)
+
+        mesh.writeVTK('temp.vtr', {'arange.txt':vec})
+        meshVTR, models = TensorMesh.readVTK('temp.vtr')
+
+        assert mesh.__str__() == meshVTR.__str__()
+        assert np.all(np.array(mesh.h) - np.array(meshVTR.h) == 0)
+
+        assert 'arange.txt' in models
+        vecVTK = models['arange.txt']
+        assert np.sum(vec - vecVTK) == 0
+
+        print 'IO of VTR tensor mesh files is working'
+        os.remove('temp.vtr')
+
+
+class TestOcTreeMeshIO(unittest.TestCase):
+
+    def setUp(self):
+        h = np.ones(16)
+        mesh = TreeMesh([h,2*h,3*h])
+        mesh.refine(3)
+        mesh._refineCell([0,0,0,3])
+        mesh._refineCell([0,2,0,3])
+        self.mesh = mesh
+
+    def test_UBCfiles(self):
+
+        mesh = self.mesh
+        # Make a vector
+        vec = np.arange(mesh.nC)
+        # Write and read
+        mesh.writeUBC('temp.msh', {'arange.txt':vec})
+        meshUBC = TreeMesh.readUBC('temp.msh')
+        vecUBC  = meshUBC.readModelUBC('arange.txt')
+
+        # The mesh
+        assert mesh.__str__() == meshUBC.__str__()
+        assert np.sum(mesh.gridCC - meshUBC.gridCC) == 0
+        assert np.sum(vec - vecUBC) == 0
+        assert np.all(np.array(mesh.h) - np.array(meshUBC.h) == 0)
+        print 'IO of UBC octree files is working'
+        os.remove('temp.msh')
+        os.remove('arange.txt')
+
+    def test_VTUfiles(self):
+        mesh = self.mesh
+        vec = np.arange(mesh.nC)
+        mesh.writeVTK('temp.vtu',{'arange':vec})
+        print 'Writing of VTU files is working'
+        os.remove('temp.vtu')
+
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/mesh/test_TreeMesh.py

@@ -26,6 +26,27 @@ class TestSimpleQuadTree(unittest.TestCase):
 
         assert np.allclose(np.r_[M._areaFxFull, M._areaFyFull], M._deflationMatrix('F') * M.area)
 
+    def test_getitem(self):
+        M = Mesh.TreeMesh([4,4])
+        M.refine(1)
+        assert M.nC == 4
+        assert len(M) == M.nC
+        assert np.allclose(M[0].center, [0.25,0.25])
+        actual = [[0,0],[0.5,0],[0,0.5],[0.5,0.5]]
+        for i, n in enumerate(M[0].nodes):
+            assert np.allclose(M._gridN[n,:], actual[i])
+
+    def test_getitem3D(self):
+        M = Mesh.TreeMesh([4,4,4])
+        M.refine(1)
+        assert M.nC == 8
+        assert len(M) == M.nC
+        assert np.allclose(M[0].center, [0.25,0.25,0.25])
+        actual = [[0,0,0],[0.5,0,0],[0,0.5,0],[0.5,0.5,0],
+                  [0,0,0.5],[0.5,0,0.5],[0,0.5,0.5],[0.5,0.5,0.5]]
+        for i, n in enumerate(M[0].nodes):
+            assert np.allclose(M._gridN[n,:], actual[i])
+
     def test_refine(self):
         M = Mesh.TreeMesh([4,4,4])
         M.refine(1)

tests/utils/test_Zero.py

@@ -1,5 +1,5 @@
 import unittest
-from SimPEG.Utils import Zero, Identity, sdiag
+from SimPEG.Utils import Zero, Identity, sdiag, mkvc
 from SimPEG import np, sp
 
 class Tests(unittest.TestCase):
@@ -29,6 +29,11 @@ class Tests(unittest.TestCase):
         assert a == 1
         self.assertRaises(ZeroDivisionError, lambda:3/z)
 
+        assert mkvc(z) == 0
+        assert sdiag(z)*a == 0
+        assert z.T == 0
+        assert z.transpose == 0
+
     def test_mat_zero(self):
         z = Zero()
         S = sdiag(np.r_[2,3])