From 83cb5ce46a9dd3e908e95935856f7910dc8bd65a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 28 Nov 2015 12:55:24 -0800 Subject: [PATCH 01/30] - Meshlesses Identity Map (takes nP instead of a mesh) - Tikhonov regularization if active cells are used (don't take derivs across interfaces between active cells and not) - testing improvements: test 1D, 2D, 3D on a random tensor mesh , also test that for a constant mref, phi_m(ref) = 0 --- SimPEG/Maps.py | 72 +++++++++++++++++++++++++++++ SimPEG/Regularization.py | 49 +++++++++++++++++--- tests/base/test_regularization.py | 75 ++++++++++++++++++++++++++----- 3 files changed, 179 insertions(+), 17 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 5b4782ac..7c46fa1d 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -118,6 +118,78 @@ class IdentityMap(object): def __str__(self): return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1]) + +class IdentityMap_Meshless(IdentityMap): + + def __init__(self, nP=None, **kwargs): + IdentityMap.__init__(self, None, **kwargs) + self._nP = nP + + @property + def nP(self): + """ + :rtype: int + :return: number of parameters in the model + """ + if self._nP is None: + return '*' + return self._nP + + @property + def shape(self): + """ + The default shape is (mesh.nC, nP). + + :rtype: (int,int) + :return: shape of the operator as a tuple + """ + if self._nP is None: + return ('*', '*') + return (self.nP, self.nP) + + + def _transform(self, m): + """ + Changes the model into the physical property. + + .. note:: + + This can be called by the __mul__ property against a numpy.ndarray. + + :param numpy.array m: model + :rtype: numpy.array + :return: transformed model + + """ + return m + + def inverse(self, D): + """ + Changes the physical property into the model. + + .. note:: + + The *transformInverse* may not be easy to create in general. + + :param numpy.array D: physical property + :rtype: numpy.array + :return: model + + """ + raise NotImplementedError('The transformInverse is not implemented.') + + def deriv(self, m): + """ + The derivative of the transformation. + + :param numpy.array m: model + :rtype: scipy.csr_matrix + :return: derivative of transformed model + + """ + return sp.identity(self.nP) + + class ComboMap(IdentityMap): """Combination of various maps.""" diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index e3506290..7e98f073 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -20,12 +20,13 @@ class BaseRegularization(object): mesh = None #: A SimPEG.Mesh instance. mref = None #: Reference model. - def __init__(self, mesh, mapping=None, **kwargs): + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Utils.setKwargs(self, **kwargs) self.mesh = mesh assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object." self.mapping = mapping or Maps.IdentityMap(mesh) self.mapping._assertMatchesPair(self.mapPair) + self.indActive = indActive @property def parent(self): @@ -112,8 +113,6 @@ class BaseRegularization(object): return mD.T * ( self.W.T * ( self.W * ( mD * v) ) ) - - class Tikhonov(BaseRegularization): """**Tikhonov Regularization** @@ -205,14 +204,18 @@ class Tikhonov(BaseRegularization): alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction") alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") - def __init__(self, mesh, mapping=None, **kwargs): + def __init__(self, mesh, mapping=None, indActive = None, **kwargs): BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs) + self.indActive = indActive @property def Ws(self): """Regularization matrix Ws""" if getattr(self,'_Ws', None) is None: - self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) + self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Ws = Pac.T * self._Ws * Pac return self._Ws @property @@ -221,6 +224,13 @@ class Tikhonov(BaseRegularization): if getattr(self, '_Wx', None) is None: Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx + + if self.indActive is not None: + indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1 + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx] + self._Wx = Pafx.T*self._Wx*Pac + return self._Wx @property @@ -229,6 +239,13 @@ class Tikhonov(BaseRegularization): if getattr(self, '_Wy', None) is None: Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady + + if self.indActive is not None: + indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1 + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy] + self._Wy = Pafy.T*self._Wy*Pac + return self._Wy @property @@ -237,6 +254,13 @@ class Tikhonov(BaseRegularization): if getattr(self, '_Wz', None) is None: Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz + + if self.indActive is not None: + indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1 + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz] + self._Wz = Pafz.T*self._Wz*Pac + return self._Wz @property @@ -244,6 +268,11 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wxx""" if getattr(self, '_Wxx', None) is None: self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx + + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Wxx = Pac.T*self._Wxx*Pac + return self._Wxx @property @@ -251,6 +280,11 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wyy""" if getattr(self, '_Wyy', None) is None: self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady + + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Wyy = Pac.T*self._Wyy*Pac + return self._Wyy @property @@ -258,6 +292,11 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wzz""" if getattr(self, '_Wzz', None) is None: self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz + + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Wzz = Pac.T*self._Wzz*Pac + return self._Wzz @property diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py index af7da692..3243e51d 100644 --- a/tests/base/test_regularization.py +++ b/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) + + for i, mesh in enumerate(self.meshlist): - 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) + 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_Meshless(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__': From cfc921b66715c6935813fece4d631d6516afee3b Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 28 Nov 2015 13:12:44 -0800 Subject: [PATCH 02/30] cleaned out transform, inverse and deriv (all are inherited from IdentityMap) --- SimPEG/Maps.py | 42 ------------------------------------------ 1 file changed, 42 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 7c46fa1d..d36deed4 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -148,48 +148,6 @@ class IdentityMap_Meshless(IdentityMap): return (self.nP, self.nP) - def _transform(self, m): - """ - Changes the model into the physical property. - - .. note:: - - This can be called by the __mul__ property against a numpy.ndarray. - - :param numpy.array m: model - :rtype: numpy.array - :return: transformed model - - """ - return m - - def inverse(self, D): - """ - Changes the physical property into the model. - - .. note:: - - The *transformInverse* may not be easy to create in general. - - :param numpy.array D: physical property - :rtype: numpy.array - :return: model - - """ - raise NotImplementedError('The transformInverse is not implemented.') - - def deriv(self, m): - """ - The derivative of the transformation. - - :param numpy.array m: model - :rtype: scipy.csr_matrix - :return: derivative of transformed model - - """ - return sp.identity(self.nP) - - class ComboMap(IdentityMap): """Combination of various maps.""" From c4d34c4e0d8083641f79a745196a51aebfb8b29d Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 30 Nov 2015 17:46:57 -0800 Subject: [PATCH 03/30] Initial counting of nodes. - Some of the nodes in the cell may be hanging. --- SimPEG/Mesh/TreeMesh.py | 62 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 59 insertions(+), 3 deletions(-) diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 1e6c91c0..9554720f 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -824,8 +824,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts): def _numberCells(self, force=False): if not self.__dirtyCells__ and not force: return self._cc2i = dict() + self._i2cc = dict() for ii, c in enumerate(sorted(self._cells)): self._cc2i[c] = ii + self._i2cc[ii] = c self.__dirtyCells__ = False def _numberNodes(self, force=False): @@ -2181,6 +2183,25 @@ class TreeMesh(BaseTensorMesh, InnerProducts): if showIt: plt.show() return tuple(out) + def __len__(self): return self.nC + + def __getitem__(self, key): + if isinstance( key, slice ) : + #Get the start, stop, and step from the slice + return [self[ii] for ii in xrange(*key.indices(len(self)))] + elif isinstance( key, int ) : + if key < 0 : #Handle negative indices + key += len( self ) + if key >= len( self ) : + raise IndexError, "The index (%d) is out of range."%key + + self._numberCells() # no-op if numbered + index = self._i2cc[key] + pointer = self._asPointer(index) + return Cell(self, index, pointer) + else: + raise TypeError, "Invalid argument type." + class Cell(object): def __init__(self, mesh, index, pointer): @@ -2188,6 +2209,34 @@ class Cell(object): self._index = index self._pointer = pointer + @property + def nodes(self): + M = self.mesh + M._numberNodes() + p = self._pointer + i = self._index + w = M._levelWidth(p[-1]) + + if M.dim == 2: + n = [ + i, + M._index([ p[0] + w, p[1] , p[2]]), + M._index([ p[0] , p[1]+ w, p[2]]), + M._index([ p[0] + w, p[1]+ w, p[2]]), + ] + elif self.dim == 3: + n = [ + i, + M._index([ p[0] + w, p[1] , p[2] ,p[3]]), + M._index([ p[0] , p[1] + w, p[2] ,p[3]]), + M._index([ p[0] + w, p[1] + w, p[2] ,p[3]]), + M._index([ p[0] , p[1] , p[2] + w,p[3]]), + M._index([ p[0] + w, p[1] , p[2] + w,p[3]]), + M._index([ p[0] , p[1] + w, p[2] + w,p[3]]), + M._index([ p[0] + w, p[1] + w, p[2] + w,p[3]]), + ] + return [M._n2i[_] for _ in n] + @property def center(self): if getattr(self, '_center', None) is None: @@ -2270,7 +2319,7 @@ if __name__ == '__main__': # T = TreeMesh([[(1,128)],[(1,128)],[(1,128)]],levels=7) # T = TreeMesh([128,128,128]) # T = TreeMesh([64,64],levels=6) - T = TreeMesh([4,4,4]) + T = TreeMesh([8,8]) # T = TreeMesh([[(1,128)],[(1,128)]],levels=7) # T.refine(lambda xc:2, balance=False) # T._index([0,0,0]) @@ -2281,8 +2330,15 @@ if __name__ == '__main__': T.refine(function)#, balance=False) # print time.time() - tic # print T.nC - T.plotSlice(np.log(T.vol))#np.random.rand(T.nC)) - + # T.plotSlice(np.log(T.vol))#np.random.rand(T.nC)) + T.plotGrid() + # print [c for c in T] + c = T[0] + plt.plot(c.center[0],c.center[1],'r.') + nodes = c.nodes + for n in nodes: + _ = T._gridN[n,:] + plt.plot(_[0],_[1],'gs') plt.show() blah From e30a7bcafc963c3ed6ebb168df3c97c70e5ce3e8 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 30 Nov 2015 17:52:38 -0800 Subject: [PATCH 04/30] documentation updates --- SimPEG/Mesh/TreeMesh.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 9554720f..a4515dab 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -561,15 +561,18 @@ class TreeMesh(BaseTensorMesh, InnerProducts): return [p - (p % mod) for p in pointer[:-1]] + [pointer[-1]-1] def _cellN(self, p): + """Node location [x,y(,z)] of a single cell, closest to origin, given a pointer.""" p = self._asPointer(p) return [hi[:p[ii]].sum() for ii, hi in enumerate(self.h)] def _cellH(self, p): + """Widths of a single cell given a pointer.""" p = self._asPointer(p) w = self._levelWidth(p[-1]) return [hi[p[ii]:p[ii]+w].sum() for ii, hi in enumerate(self.h)] def _cellC(self, p): + """Cell center of a single cell (without origin correction), given a pointer.""" return (np.array(self._cellH(p))/2.0 + self._cellN(p)).tolist() def _levelWidth(self, level): From a7ab0dc1e246eabb440a56a3e25bbd5851e950ad Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 30 Nov 2015 18:00:16 -0800 Subject: [PATCH 05/30] Unit tests for getitem on tree mesh --- tests/mesh/test_TreeMesh.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/mesh/test_TreeMesh.py b/tests/mesh/test_TreeMesh.py index e624ce87..afad27d1 100644 --- a/tests/mesh/test_TreeMesh.py +++ b/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) From 7da637e883304ca1ecbaf016357a0b3d51bf9a1b Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 30 Nov 2015 18:05:02 -0800 Subject: [PATCH 06/30] documentation on Cell.nodes --- SimPEG/Mesh/TreeMesh.py | 1 + 1 file changed, 1 insertion(+) diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index a4515dab..eea2ad4c 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -2214,6 +2214,7 @@ class Cell(object): @property def nodes(self): + """The node index in _gridN (this may include hanging nodes).""" M = self.mesh M._numberNodes() p = self._pointer From 589cd655af735394e7dc322ad496eb04179b8179 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 2 Dec 2015 16:00:01 -0800 Subject: [PATCH 07/30] Updated vtk write classes. --- SimPEG/Utils/meshutils.py | 80 +++++++++++++++++++++++++++++++++++---- 1 file changed, 72 insertions(+), 8 deletions(-) diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 585fcc9a..816132fc 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -205,7 +205,6 @@ def writeUBCTensorModel(fileName, mesh, model): np.savetxt(fileName, modelMatTR.ravel()) - def readVTRFile(fileName): """ Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model @@ -296,13 +295,14 @@ def writeVTRFile(fileName,mesh,model=None): vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1)) # Assign the model('s) to the object - for item in model.iteritems(): - # Convert numpy array - vtkDoubleArr = numpy_to_vtk(item[1],deep=1) - vtkDoubleArr.SetName(item[0]) - vtkObj.GetCellData().AddArray(vtkDoubleArr) - # Set the active scalar - vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) + if model is not None: + for item in model.iteritems(): + # Convert numpy array + vtkDoubleArr = numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtkObj.GetCellData().AddArray(vtkDoubleArr) + # Set the active scalar + vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) vtkObj.Update() @@ -318,6 +318,70 @@ def writeVTRFile(fileName,mesh,model=None): vtrWriteFilter.SetFileName(fileName) vtrWriteFilter.Update() +def writeVTUFile(fileName,ocTreeMesh,modelDict=None): + ''' + Function to write a VTU file from a SimPEG TreeMesh and model. + ''' + from vtk import vtkXMLUnstructuredGridWriter as Writer + from vtk.util.numpy_support import numpy_to_vtk + + # Make the object + vtuObj = simpegOcTree2vtuObj(ocTreeMesh,modelDict) + + # Make the writer + vtuWriteFilter = Writer() + if float(vtk.VTK_VERSION.split('.')[0]) >=6: + vtuWriteFilter.SetInputData(vtuObj) + else: + vtuWriteFilter.SetInput(vtuObj) + vtuWriteFilter.SetInput(vtuObj) + vtuWriteFilter.SetFileName(fileName) + # Write the file + vtuWriteFilter.Update() + +def simpegOcTree2vtuObj(simpegOcTreeMesh,modelDict=None): + ''' + Convert simpeg OcTree mesh and model to a VTK vtu object. + + ''' + import vtk + from vtk.util.numpy_support import numpy_to_vtk + + if str(type(simpegOcTreeMesh)).split()[-1][1:-2] not in 'SimPEG.Mesh.TreeMesh.TreeMesh': + raise IOError('simpegOcTreeMesh is not a SimPEG TreeMesh.') + + # Make the data parts for the vtu object + # Points + ptsMat = simpegOcTreeMesh._gridN + simpegOcTreeMesh.x0 + vtkPts = vtk.vtkPoints() + vtkPts.SetData(npsup.numpy_to_vtk(ptsMat,deep=True)) + # Cells + cellConn = np.array([c.nodes for c in simpegOcTreeMesh],dtype=np.int64) + + cellsMat = np.concatenate((np.ones((cellConn.shape[0],1),dtype=np.int64)*cellConn.shape[1],cellConn),axis=1).ravel() + cellsArr = vtk.vtkCellArray() + cellsArr.SetNumberOfCells(cellConn.shape[0]) + cellsArr.SetCells(cellConn.shape[0],npsup.numpy_to_vtkIdTypeArray(cellsMat,deep=True)) + + # Make the object + vtuObj = vtk.vtkUnstructuredGrid() + vtuObj.SetPoints(vtkPts) + vtuObj.SetCells(vtk.VTK_VOXEL,cellsArr) + # Add the level of refinement as a cell array + cellSides = np.array([np.array(vtuObj.GetCell(i).GetBounds()).reshape((3,2)).dot(np.array([-1, 1])) for i in np.arange(vtuObj.GetNumberOfCells())]) + uniqueLevel, indLevel = np.unique(np.prod(cellSides,axis=1),return_inverse=True) + refineLevelArr = npsup.numpy_to_vtk(indLevel.max() - indLevel,deep=1) + refineLevelArr.SetName('octreeLevel') + vtuObj.GetCellData().AddArray(refineLevelArr) + # Assign the model('s) to the object + if modelDict is not None: + for item in modelDict.iteritems(): + # Convert numpy array + vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtuObj.GetCellData().AddArray(vtkDoubleArr) + + return vtuObj def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'): """ From a8551f3e04473a94e4e2c3bf044a0f515859c9bb Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 2 Dec 2015 16:04:28 -0800 Subject: [PATCH 08/30] Added function to write a UBC octree mesh for TreeMesh object. --- SimPEG/Utils/meshutils.py | 51 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 816132fc..81b15566 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -205,6 +205,57 @@ def writeUBCTensorModel(fileName, mesh, model): np.savetxt(fileName, modelMatTR.ravel()) +def writeUBCocTreeFiles(fileName,mesh,modelDict=None,): + ''' + Write UBC ocTree mesh and model files from a simpeg ocTree mesh and model. + + :param + + ''' + + # Calculate information to write in the file. + # Number of cells in the underlying mesh + nCunderMesh = np.array([h.size for h in mesh.h],dtype=np.int64) + # The top-south-west most corner of the mesh + tswCorn = mesh.x0 + np.array([0,0,np.sum(mesh.h[2])]) + # Smallest cell size + smallCell = np.array([h.min() for h in mesh.h]) + # Number of cells + nrCells = mesh.nC + + ## Extract iformation about the cells. + # cell pointers + cellPointers = np.array([c._pointer for c in mesh]) + # cell with + cellW = np.array([ mesh._levelWidth(i) for i in cellPointers[:,-1] ]) + # Need to shift the pointers to work with UBC indexing + # UBC Octree indexes always the top-left-close (top-south-west) corner first and orders the cells in z(top-down),x,y vs x,y,z(bottom-up). + # Shift index up by 1 + ubcCellPt = cellPointers[:,0:-1].copy() + np.array([1.,1.,1.]) + # Need reindex the z index to be from the top-left-close corner and to be from the global top. + ubcCellPt[:,2] = ( nCunderMesh[-1] + 2) - (ubcCellPt[:,2] + cellW) + + # Reorder the ubcCellPt + ubcReorder = np.argsort(ubcCellPt.view(','.join(3*['float'])),axis=0,order=['f2','f1','f0'])[:,0] + # Make a array with the pointers and the withs, that are order in the ubc ordering + indArr = np.concatenate((ubcCellPt[ubcReorder,:],cellW[ubcReorder].reshape((-1,1)) ),axis=1) + + ## Write the UBC octree mesh file + with open(fileName,'w') as mshOut: + mshOut.write('{:.0f} {:.0f} {:.0f}\n'.format(nCunderMesh[0],nCunderMesh[1],nCunderMesh[2])) + mshOut.write('{:.4f} {:.4f} {:.4f}\n'.format(tswCorn[0],tswCorn[1],tswCorn[2])) + mshOut.write('{:.3f} {:.3f} {:.3f}\n'.format(smallCell[0],smallCell[1],smallCell[2])) + mshOut.write('{:.0f} \n'.format(nrCells)) + np.savetxt(mshOut,indArr,fmt='%i') + + ## Print the models + # Assign the model('s) to the object + if modelDict is not None: + # indUBCvector = np.argsort(cX0[np.argsort(np.concatenate((cX0[:,0:2],cX0[:,2:3].max() - cX0[:,2:3]),axis=1).view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0]].view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0] + for item in modelDict.iteritems(): + # Save the data + np.savetxt(item[0],item[1][ubcReorder],fmt='%3.5e') + def readVTRFile(fileName): """ Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model From 25ad1488f5c617a2f9d24aab1af93db76d372595 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 2 Dec 2015 19:20:21 -0800 Subject: [PATCH 09/30] Added a function to read UBC octree mesh. Updated __init__ to import the new functions. --- SimPEG/Utils/__init__.py | 2 +- SimPEG/Utils/meshutils.py | 66 +++++++++++++++++++++++++++++++++++++-- 2 files changed, 65 insertions(+), 3 deletions(-) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 5280ae79..4d5c8e2e 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, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile +from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCocTreeFiles, readUBCocTreeFiles, writeUBCTensorModel, readVTRFile, writeVTRFile, writeVTUFile from curvutils 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 81b15566..3537f0b5 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -205,11 +205,13 @@ def writeUBCTensorModel(fileName, mesh, model): np.savetxt(fileName, modelMatTR.ravel()) -def writeUBCocTreeFiles(fileName,mesh,modelDict=None,): +def writeUBCocTreeFiles(fileName,mesh,modelDict=None): ''' Write UBC ocTree mesh and model files from a simpeg ocTree mesh and model. - :param + :param str fileName: File to write to + :param simpeg.Mesh.TreeMesh mesh: The mesh + :param dictionary modelDict: The models in a dictionary, where the keys is the name of the of the model file ''' @@ -256,6 +258,66 @@ def writeUBCocTreeFiles(fileName,mesh,modelDict=None,): # Save the data np.savetxt(item[0],item[1][ubcReorder],fmt='%3.5e') +def readUBCocTreeFiles(meshFile,modelFiles=None): + """ + Read UBC 3D OcTree mesh and/or modelFiles + + Input: + :param str meshFile: path to the UBC GIF OcTree mesh file to read + :param list of str modelFiles: list of paths modelFiles + + Output: + :return SimPEG.Mesh.TreeMesh mesh: The octree mesh + :return list of ndarray's: models as a list of numpy array's + """ + + ## Read the file lines + fileLines = np.genfromtxt(meshFile,dtype=str,delimiter='\n') + # Extract the data + nCunderMesh = np.array(fileLines[0].split(),dtype=float) + # I think this is the case? + if np.unique(nCunderMesh).size >1: + raise Exception('SimPEG TreeMeshes have the same number of cell in all directions') + tswCorn = np.array(fileLines[1].split(),dtype=float) + smallCell = np.array(fileLines[2].split(),dtype=float) + nrCells = np.array(fileLines[3].split(),dtype=float) + # Read the index array + indArr = np.genfromtxt(fileLines[4::],dtype=np.int) + + ## Calculate simpeg parameters + h1,h2,h3 = [np.ones(nr)*sz for nr,sz in zip(nCunderMesh,smallCell)] + x0 = tswCorn - np.sum(h3) + # Need to convert the index array to a points list that complies with SimPEG TreeMesh. + # Shift to start at 0 + simpegCellPt = indArr[:,0:-1].copy() - np.array([1.,1.,1.]) + # Need reindex the z index to be from the bottom-left-close corner and to be from the global bottom. + simpegCellPt[:,2] = ( nCunderMesh[-1] + 2) - (simpegCellPt[:,2] - indArr[:,3]) + # Figure out the reordering + simpegReorder = np.argsort(simpegCellPt.view(','.join(3*['float'])),axis=0,order=['f2','f1','f0'])[:,0] + + # Calculate the cell level + simpegLevel = np.log2(np.min(nCunderMesh)) - np.log2(indArr[:,3]) + # Make a pointer matrix + simpegPointers = np.concatenate((simpegCellPt[simpegReorder,:],simpegLevel[simpegReorder].reshape((-1,1))),axis=1) + # Make an index set + + ## Make the tree mesh + from SimPEG.Mesh import TreeMesh + mesh = TreeMesh([h1,h2,h3],x0) + mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()]) + + if modelFiles is None: + return mesh + else: + modList = [] + for modFile in modelFiles: + modArr = np.loadtxt(modFile) + if len(modArr.shape) == 1: + modList.append(modArr[simpegReorder]) + else: + modList.append(modArr[simpegReorder,:]) + return mesh, modList + def readVTRFile(fileName): """ Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model From c298ebe8d8fd33db12d97c8d9013a1148dfd0449 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 4 Dec 2015 15:42:08 -0800 Subject: [PATCH 10/30] Remove the Meshless Identity Map. - This is now default functionality in the IdentityMap. --- SimPEG/Maps.py | 127 +++++++++++++----------------- tests/base/test_regularization.py | 8 +- 2 files changed, 57 insertions(+), 78 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index d36deed4..a3c76e6a 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -10,21 +10,25 @@ class IdentityMap(object): SimPEG Map """ - __metaclass__ = Utils.SimPEGMetaClass - mesh = None #: A SimPEG Mesh - - def __init__(self, mesh, **kwargs): + def __init__(self, mesh=None, nP=None, **kwargs): Utils.setKwargs(self, **kwargs) + + if nP is not None: + assert type(nP) in [int, long], ' Number of parameters must be an integer.' + self.mesh = mesh + self._nP = nP @property def nP(self): """ :rtype: int - :return: number of parameters in the model + :return: number of parameters that the mapping accepts """ + if self._nP is not None: + return self._nP if self.mesh is None: return '*' return self.mesh.nC @@ -32,11 +36,15 @@ class IdentityMap(object): @property def shape(self): """ - The default shape is (mesh.nC, nP). + The default shape is (mesh.nC, nP) if the mesh is defined. + If this is a meshless mapping (i.e. nP is defined independently) + the shape will be the the shape (nP,nP). :rtype: (int,int) :return: shape of the operator as a tuple """ + if self._nP is not None: + return (self.nP, self.nP) if self.mesh is None: return ('*', self.nP) return (self.mesh.nC, self.nP) @@ -119,35 +127,6 @@ class IdentityMap(object): return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1]) -class IdentityMap_Meshless(IdentityMap): - - def __init__(self, nP=None, **kwargs): - IdentityMap.__init__(self, None, **kwargs) - self._nP = nP - - @property - def nP(self): - """ - :rtype: int - :return: number of parameters in the model - """ - if self._nP is None: - return '*' - return self._nP - - @property - def shape(self): - """ - The default shape is (mesh.nC, nP). - - :rtype: (int,int) - :return: shape of the operator as a tuple - """ - if self._nP is None: - return ('*', '*') - return (self.nP, self.nP) - - class ComboMap(IdentityMap): """Combination of various maps.""" @@ -505,7 +484,7 @@ class ActiveCells(IdentityMap): else: self.valInactive = valInactive.copy() self.valInactive[self.indActive] = 0 - + inds = np.nonzero(self.indActive)[0] self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP)) @@ -738,7 +717,7 @@ class PolyMap(IdentityMap): Parameterize the model space using a polynomials in a wholespace. ..math:: - + y = \mathbf{V} c Define the model as: @@ -782,10 +761,10 @@ class PolyMap(IdentityMap): else: raise(Exception("Input for normal = X or Y or Z")) #3D - elif self.mesh.dim == 3: + elif self.mesh.dim == 3: X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] - Z = self.mesh.gridCC[:,2] + Y = self.mesh.gridCC[:,1] + Z = self.mesh.gridCC[:,2] if self.normal =='X': f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X elif self.normal =='Y': @@ -796,43 +775,43 @@ class PolyMap(IdentityMap): raise(Exception("Input for normal = X or Y or Z")) else: raise(Exception("Only supports 2D")) - + return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5) - + def deriv(self, m): alpha = self.slope sig1,sig2, c = m[0],m[1],m[2:] if self.logSigma: sig1, sig2 = np.exp(sig1), np.exp(sig2) #2D - if self.mesh.dim == 2: + if self.mesh.dim == 2: X = self.mesh.gridCC[:,0] Y = self.mesh.gridCC[:,1] if self.normal =='X': f = polynomial.polyval(Y, c) - X - V = polynomial.polyvander(Y, len(c)-1) + V = polynomial.polyvander(Y, len(c)-1) elif self.normal =='Y': f = polynomial.polyval(X, c) - Y - V = polynomial.polyvander(X, len(c)-1) + V = polynomial.polyvander(X, len(c)-1) else: - raise(Exception("Input for normal = X or Y or Z")) + raise(Exception("Input for normal = X or Y or Z")) #3D - elif self.mesh.dim == 3: + elif self.mesh.dim == 3: X = self.mesh.gridCC[:,0] Y = self.mesh.gridCC[:,1] Z = self.mesh.gridCC[:,2] if self.normal =='X': f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X - V = polynomial.polyvander2d(Y, Z, self.order) + V = polynomial.polyvander2d(Y, Z, self.order) elif self.normal =='Y': f = polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - Y - V = polynomial.polyvander2d(X, Z, self.order) + V = polynomial.polyvander2d(X, Z, self.order) elif self.normal =='Z': f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z - V = polynomial.polyvander2d(X, Y, self.order) + V = polynomial.polyvander2d(X, Y, self.order) else: raise(Exception("Input for normal = X or Y or Z")) @@ -845,16 +824,16 @@ class PolyMap(IdentityMap): g3 = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*V - return sp.csr_matrix(np.c_[g1,g2,g3]) + return sp.csr_matrix(np.c_[g1,g2,g3]) class SplineMap(IdentityMap): """SplineMap - Parameterize the boundary of two geological units using a spline interpolation + Parameterize the boundary of two geological units using a spline interpolation ..math:: - + g = f(x)-y Define the model as: @@ -879,7 +858,7 @@ class SplineMap(IdentityMap): def nP(self): if self.mesh.dim == 2: return np.size(self.pts)+2 - elif self.mesh.dim == 3: + elif self.mesh.dim == 3: return np.size(self.pts)*2+2 else: raise(Exception("Only supports 2D and 3D")) @@ -896,28 +875,28 @@ class SplineMap(IdentityMap): X = self.mesh.gridCC[:,0] Y = self.mesh.gridCC[:,1] self.spl = UnivariateSpline(self.pts, c, k=self.order, s=0) - if self.normal =='X': + if self.normal =='X': f = self.spl(Y) - X elif self.normal =='Y': f = self.spl(X) - Y else: raise(Exception("Input for normal = X or Y or Z")) - # 3D: - # Comments: + # 3D: + # Comments: # Make two spline functions and link them using linear interpolation. # This is not quite direct extension of 2D to 3D case # Using 2D interpolation is possible - elif self.mesh.dim == 3: + elif self.mesh.dim == 3: X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] + Y = self.mesh.gridCC[:,1] Z = self.mesh.gridCC[:,2] - npts = np.size(self.pts) + npts = np.size(self.pts) if np.mod(c.size, 2): raise(Exception("Put even points!")) - + self.spl = {"splb":UnivariateSpline(self.pts, c[:npts], k=self.order, s=0), "splt":UnivariateSpline(self.pts, c[npts:], k=self.order, s=0)} @@ -932,7 +911,7 @@ class SplineMap(IdentityMap): raise(Exception("Input for normal = X or Y or Z")) else: raise(Exception("Only supports 2D and 3D")) - + return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5) @@ -942,7 +921,7 @@ class SplineMap(IdentityMap): if self.logSigma: sig1, sig2 = np.exp(sig1), np.exp(sig2) #2D - if self.mesh.dim == 2: + if self.mesh.dim == 2: X = self.mesh.gridCC[:,0] Y = self.mesh.gridCC[:,1] @@ -951,9 +930,9 @@ class SplineMap(IdentityMap): elif self.normal =='Y': f = self.spl(X) - Y else: - raise(Exception("Input for normal = X or Y or Z")) + raise(Exception("Input for normal = X or Y or Z")) #3D - elif self.mesh.dim == 3: + elif self.mesh.dim == 3: X = self.mesh.gridCC[:,0] Y = self.mesh.gridCC[:,1] Z = self.mesh.gridCC[:,2] @@ -961,7 +940,7 @@ class SplineMap(IdentityMap): zb = self.ptsv[0] zt = self.ptsv[1] flines = (self.spl["splt"](Y)-self.spl["splb"](Y))*(Z-zb)/(zt-zb) + self.spl["splb"](Y) - f = flines - X + f = flines - X # elif self.normal =='Y': # elif self.normal =='Z': else: @@ -974,7 +953,7 @@ class SplineMap(IdentityMap): g1 = -(np.arctan(alpha*f)/np.pi + 0.5) + 1.0 g2 = (np.arctan(alpha*f)/np.pi + 0.5) - + if self.mesh.dim ==2: g3 = np.zeros((self.mesh.nC, self.npts)) if self.normal =='Y': @@ -988,7 +967,7 @@ class SplineMap(IdentityMap): cb = c.copy() dy = self.mesh.hy[ind]*1.5 ca[i] = ctemp+dy - cb[i] = ctemp-dy + cb[i] = ctemp-dy spla = UnivariateSpline(self.pts, ca, k=self.order, s=0) splb = UnivariateSpline(self.pts, cb, k=self.order, s=0) fderiv = (spla(X)-splb(X))/(2*dy) @@ -998,7 +977,7 @@ class SplineMap(IdentityMap): g3 = np.zeros((self.mesh.nC, self.npts*2)) if self.normal =='X': # Here we use perturbation to compute sensitivity - for i in range(self.npts*2): + for i in range(self.npts*2): ctemp = c[i] ind = np.argmin(abs(self.mesh.vectorCCy-ctemp)) ca = c.copy() @@ -1012,20 +991,20 @@ class SplineMap(IdentityMap): splbb = UnivariateSpline(self.pts, cb[:self.npts], k=self.order, s=0) flinesa = (self.spl["splt"](Y)-splba(Y))*(Z-zb)/(zt-zb) + splba(Y) - X flinesb = (self.spl["splt"](Y)-splbb(Y))*(Z-zb)/(zt-zb) + splbb(Y) - X - #treat top boundary + #treat top boundary else: splta = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0) spltb = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0) flinesa = (self.spl["splt"](Y)-splta(Y))*(Z-zb)/(zt-zb) + splta(Y) - X - flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X - fderiv = (flinesa-flinesb)/(2*dy) + flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X + fderiv = (flinesa-flinesb)/(2*dy) g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv else : raise(Exception("Not Implemented for Y and Z, your turn :)")) - return sp.csr_matrix(np.c_[g1,g2,g3]) + return sp.csr_matrix(np.c_[g1,g2,g3]) + - \ No newline at end of file diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py index 3243e51d..050c46ac 100644 --- a/tests/base/test_regularization.py +++ b/tests/base/test_regularization.py @@ -22,7 +22,7 @@ class RegularizationTests(unittest.TestCase): if not inspect.isclass(r): continue if not issubclass(r, Regularization.BaseRegularization): continue - + for i, mesh in enumerate(self.meshlist): print 'Testing %iD'%mesh.dim @@ -32,7 +32,7 @@ class RegularizationTests(unittest.TestCase): 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) + print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) passed = reg.eval(reg.mref) < TOL self.assertTrue(passed) @@ -50,7 +50,7 @@ class RegularizationTests(unittest.TestCase): 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) @@ -62,7 +62,7 @@ class RegularizationTests(unittest.TestCase): 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_Meshless(nP=indAct.nonzero()[0].size) + mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size) reg = r(mesh, mapping=mapping, indActive=indAct) m = np.random.rand(mesh.nC)[indAct] From 1a40e35c269e1a4974507b9141a85bece85bb5df Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 9 Dec 2015 14:39:00 -0800 Subject: [PATCH 11/30] Fixed import bugs. --- SimPEG/Utils/meshutils.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 3537f0b5..2750d8bf 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -458,7 +458,7 @@ def simpegOcTree2vtuObj(simpegOcTreeMesh,modelDict=None): ''' import vtk - from vtk.util.numpy_support import numpy_to_vtk + from vtk.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray if str(type(simpegOcTreeMesh)).split()[-1][1:-2] not in 'SimPEG.Mesh.TreeMesh.TreeMesh': raise IOError('simpegOcTreeMesh is not a SimPEG TreeMesh.') @@ -467,14 +467,14 @@ def simpegOcTree2vtuObj(simpegOcTreeMesh,modelDict=None): # Points ptsMat = simpegOcTreeMesh._gridN + simpegOcTreeMesh.x0 vtkPts = vtk.vtkPoints() - vtkPts.SetData(npsup.numpy_to_vtk(ptsMat,deep=True)) + vtkPts.SetData(numpy_to_vtk(ptsMat,deep=True)) # Cells cellConn = np.array([c.nodes for c in simpegOcTreeMesh],dtype=np.int64) cellsMat = np.concatenate((np.ones((cellConn.shape[0],1),dtype=np.int64)*cellConn.shape[1],cellConn),axis=1).ravel() cellsArr = vtk.vtkCellArray() cellsArr.SetNumberOfCells(cellConn.shape[0]) - cellsArr.SetCells(cellConn.shape[0],npsup.numpy_to_vtkIdTypeArray(cellsMat,deep=True)) + cellsArr.SetCells(cellConn.shape[0],numpy_to_vtkIdTypeArray(cellsMat,deep=True)) # Make the object vtuObj = vtk.vtkUnstructuredGrid() @@ -483,14 +483,14 @@ def simpegOcTree2vtuObj(simpegOcTreeMesh,modelDict=None): # Add the level of refinement as a cell array cellSides = np.array([np.array(vtuObj.GetCell(i).GetBounds()).reshape((3,2)).dot(np.array([-1, 1])) for i in np.arange(vtuObj.GetNumberOfCells())]) uniqueLevel, indLevel = np.unique(np.prod(cellSides,axis=1),return_inverse=True) - refineLevelArr = npsup.numpy_to_vtk(indLevel.max() - indLevel,deep=1) + refineLevelArr = numpy_to_vtk(indLevel.max() - indLevel,deep=1) refineLevelArr.SetName('octreeLevel') vtuObj.GetCellData().AddArray(refineLevelArr) # Assign the model('s) to the object if modelDict is not None: for item in modelDict.iteritems(): # Convert numpy array - vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1) + vtkDoubleArr = numpy_to_vtk(item[1],deep=1) vtkDoubleArr.SetName(item[0]) vtuObj.GetCellData().AddArray(vtkDoubleArr) From e678affe41be5e41588045b21acde0daf50681e3 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 9 Dec 2015 15:33:05 -0800 Subject: [PATCH 12/30] Changed ave[F/E]2CC to be a csr not a css, which doesn't support indexing. --- SimPEG/Mesh/TreeMesh.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index eea2ad4c..c400b956 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -1706,9 +1706,9 @@ class TreeMesh(BaseTensorMesh, InnerProducts): "Construct the averaging operator on cell faces to cell centers." if getattr(self, '_aveF2CC', None) is None: if self.dim == 2: - self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC]) + self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC]).tocsr() elif self.dim == 3: - self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]) + self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]).tocsr() return self._aveF2CC @property @@ -1716,9 +1716,9 @@ class TreeMesh(BaseTensorMesh, InnerProducts): "Construct the averaging operator on cell faces to cell centers." if getattr(self, '_aveF2CCV', None) is None: if self.dim == 2: - self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC]) + self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC]).tocsr() elif self.dim == 3: - self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]) + self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]).tocsr() return self._aveF2CCV @property From d6585dcfcdd89ce38a5664fa4023a16774c789ce Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 9 Dec 2015 15:43:51 -0800 Subject: [PATCH 13/30] Add Saving Directive --- SimPEG/Directives.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 48d7abcf..349432e4 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -206,6 +206,32 @@ class SaveOutputEveryIteration(_SaveEveryIteration): f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f)) f.close() +class SaveOutputDictEveryIteration(_SaveEveryIteration): + """SaveOutputDictEveryIteration""" + + def initialize(self): + print "SimPEG.SaveOutputDictEveryIteration will save your inversion pro + + def endIter(self): + # Save the data. + ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.r + phi_ms = 0.5*ms.dot(ms) + if self.reg.smoothModel == True: + mref = self.reg.mref + else: + mref = 0 + mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) + phi_mx = 0.5 * mx.dot(mx) + if self.prob.mesh.dim==2: + my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mr + phi_my = 0.5 * my.dot(my) + else: + phi_my = 'NaN' + if self.prob.mesh.dim==3: + mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mr + phi_mz = 0.5 * mz.dot(mz) + else: + phi_mz = 'NaN' From 4b7f7c3c147403004e00edc27d42f42714c4c5ff Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 9 Dec 2015 15:46:57 -0800 Subject: [PATCH 14/30] Fixed a spelling error --- SimPEG/Directives.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 349432e4..9e4626ca 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -210,7 +210,8 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): """SaveOutputDictEveryIteration""" def initialize(self): - print "SimPEG.SaveOutputDictEveryIteration will save your inversion pro + print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'" + def endIter(self): # Save the data. From ba2ac747409a12fb09c202d8fd61d9aa4a4ffe99 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 9 Dec 2015 15:53:41 -0800 Subject: [PATCH 15/30] Fixed code error --- SimPEG/Directives.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 9e4626ca..e91ffb71 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -210,31 +210,34 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): """SaveOutputDictEveryIteration""" def initialize(self): - print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'" - + print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName def endIter(self): # Save the data. - ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.r + ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) ) phi_ms = 0.5*ms.dot(ms) if self.reg.smoothModel == True: mref = self.reg.mref else: mref = 0 - mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) + mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) ) phi_mx = 0.5 * mx.dot(mx) if self.prob.mesh.dim==2: - my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mr + my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) ) phi_my = 0.5 * my.dot(my) else: phi_my = 'NaN' if self.prob.mesh.dim==3: - mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mr + mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) ) phi_mz = 0.5 * mz.dot(mz) else: phi_mz = 'NaN' + # Save the file as a npz + np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel) + + # class UpdateReferenceModel(Parameter): From 69d109524e50c27e29b61d02e29adbd5f8c4fdbc Mon Sep 17 00:00:00 2001 From: GudniRos Date: Sat, 12 Dec 2015 13:52:26 -0800 Subject: [PATCH 16/30] Working on reordering of UBC models. --- .gitignore | 1 + SimPEG/Utils/meshutils.py | 33 ++++++++++++++++++++------------- tests/mesh/test_MeshIO.py | 0 3 files changed, 21 insertions(+), 13 deletions(-) create mode 100644 tests/mesh/test_MeshIO.py diff --git a/.gitignore b/.gitignore index fc798cea..6a826230 100644 --- a/.gitignore +++ b/.gitignore @@ -39,3 +39,4 @@ nosetests.xml *.sublime-workspace docs/_build/ Makefile +diff.temp diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 2750d8bf..ce7e1228 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -286,25 +286,26 @@ def readUBCocTreeFiles(meshFile,modelFiles=None): ## Calculate simpeg parameters h1,h2,h3 = [np.ones(nr)*sz for nr,sz in zip(nCunderMesh,smallCell)] - x0 = tswCorn - np.sum(h3) + x0 = tswCorn - np.array([0,0,np.sum(h3)]) # Need to convert the index array to a points list that complies with SimPEG TreeMesh. # Shift to start at 0 - simpegCellPt = indArr[:,0:-1].copy() - np.array([1.,1.,1.]) + simpegCellPt = indArr[:,0:-1].copy() + simpegCellPt[:,2] = ( nCunderMesh[-1] + 2) - (simpegCellPt[:,2] + indArr[:,3]) # Need reindex the z index to be from the bottom-left-close corner and to be from the global bottom. - simpegCellPt[:,2] = ( nCunderMesh[-1] + 2) - (simpegCellPt[:,2] - indArr[:,3]) - # Figure out the reordering - simpegReorder = np.argsort(simpegCellPt.view(','.join(3*['float'])),axis=0,order=['f2','f1','f0'])[:,0] + simpegCellPt = simpegCellPt - np.array([1.,1.,1.]) # Calculate the cell level simpegLevel = np.log2(np.min(nCunderMesh)) - np.log2(indArr[:,3]) # Make a pointer matrix - simpegPointers = np.concatenate((simpegCellPt[simpegReorder,:],simpegLevel[simpegReorder].reshape((-1,1))),axis=1) - # Make an index set + simpegPointers = np.concatenate((simpegCellPt,simpegLevel.reshape((-1,1))),axis=1) + + # Figure out the reordering + simpegReorder = np.argsort((np.array([[1,1,1,-1]])*simpegPointers).view(','.join(4*['float'])),axis=0,order=['f3','f2','f1','f0'])[:,0] ## Make the tree mesh from SimPEG.Mesh import TreeMesh mesh = TreeMesh([h1,h2,h3],x0) - mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()]) + mesh._cells = set([mesh._index(p) for p in simpegPointers[simpegReorder,:].tolist()]) if modelFiles is None: return mesh @@ -427,7 +428,10 @@ def writeVTRFile(fileName,mesh,model=None): raise IOError('{:s} is an incorrect extension, has to be .vtr') # Write the file. vtrWriteFilter = rectWriter() - vtrWriteFilter.SetInput(vtkObj) + if float(VTK_VERSION.split('.')[0]) >=6: + vtrWriteFilter.SetInputData(vtkObj) + else: + vtuWriteFilter.SetInput(vtuObj) vtrWriteFilter.SetFileName(fileName) vtrWriteFilter.Update() @@ -435,7 +439,7 @@ def writeVTUFile(fileName,ocTreeMesh,modelDict=None): ''' Function to write a VTU file from a SimPEG TreeMesh and model. ''' - from vtk import vtkXMLUnstructuredGridWriter as Writer + from vtk import vtkXMLUnstructuredGridWriter as Writer, VTK_VERSION from vtk.util.numpy_support import numpy_to_vtk # Make the object @@ -443,11 +447,10 @@ def writeVTUFile(fileName,ocTreeMesh,modelDict=None): # Make the writer vtuWriteFilter = Writer() - if float(vtk.VTK_VERSION.split('.')[0]) >=6: + if float(VTK_VERSION.split('.')[0]) >=6: vtuWriteFilter.SetInputData(vtuObj) else: vtuWriteFilter.SetInput(vtuObj) - vtuWriteFilter.SetInput(vtuObj) vtuWriteFilter.SetFileName(fileName) # Write the file vtuWriteFilter.Update() @@ -465,7 +468,11 @@ def simpegOcTree2vtuObj(simpegOcTreeMesh,modelDict=None): # Make the data parts for the vtu object # Points - ptsMat = simpegOcTreeMesh._gridN + simpegOcTreeMesh.x0 + try: + ptsMat = simpegOcTreeMesh._gridN + simpegOcTreeMesh.x0 + except: + simpegOcTreeMesh.number() + ptsMat = simpegOcTreeMesh._gridN + simpegOcTreeMesh.x0 vtkPts = vtk.vtkPoints() vtkPts.SetData(numpy_to_vtk(ptsMat,deep=True)) # Cells diff --git a/tests/mesh/test_MeshIO.py b/tests/mesh/test_MeshIO.py new file mode 100644 index 00000000..e69de29b From 84eb69f6268c1d9a6bc8b0fab36d270e1d1192a9 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Mon, 14 Dec 2015 01:46:01 -0800 Subject: [PATCH 17/30] UBC ocTree read and write working. --- SimPEG/Utils/meshutils.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index ce7e1228..c43dfd13 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -299,13 +299,14 @@ def readUBCocTreeFiles(meshFile,modelFiles=None): # Make a pointer matrix simpegPointers = np.concatenate((simpegCellPt,simpegLevel.reshape((-1,1))),axis=1) - # Figure out the reordering - simpegReorder = np.argsort((np.array([[1,1,1,-1]])*simpegPointers).view(','.join(4*['float'])),axis=0,order=['f3','f2','f1','f0'])[:,0] - ## Make the tree mesh from SimPEG.Mesh import TreeMesh mesh = TreeMesh([h1,h2,h3],x0) - mesh._cells = set([mesh._index(p) for p in simpegPointers[simpegReorder,:].tolist()]) + mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()]) + + # Figure out the reordering + simpegReorder = np.argsort(np.array([mesh._index(i) for i in simpegPointers.tolist()])) + # simpegReorder = np.argsort((np.array([[1,1,1,-1]])*simpegPointers).view(','.join(4*['float'])),axis=0,order=['f3','f2','f1','f0'])[:,0] if modelFiles is None: return mesh From e42727610aad540694977adaf0f53a374c29ac90 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Mon, 14 Dec 2015 19:08:35 -0800 Subject: [PATCH 18/30] Implemented IO test for octree mesh. --- tests/mesh/test_MeshIO.py | 47 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/tests/mesh/test_MeshIO.py b/tests/mesh/test_MeshIO.py index e69de29b..e3e49548 100644 --- a/tests/mesh/test_MeshIO.py +++ b/tests/mesh/test_MeshIO.py @@ -0,0 +1,47 @@ +import numpy as np +import unittest +import SimPEG as simpeg +from SimPEG.Mesh import TensorMesh, TreeMesh + + +class TestOcTreeIO(unittest.TestCase): + + def setUp(self): + h = np.ones(16) + mesh = simpeg.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 aand read + simpeg.Utils.meshutils.writeUBCocTreeFiles('temp.msh',mesh,{'arange.txt':vec}) + meshUBC, vecUBC = simpeg.Utils.meshutils.readUBCocTreeFiles('temp.msh',['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' + + def test_VTUfiles(self): + mesh = self.mesh + vec = np.arange(mesh.nC) + try: + simpeg.Utils.meshutils.writeVTUFile('test.vtu',mesh,{'arange':vec}) + run = True + except: + run = False + assert run + print 'Writing of VTU files is working' + + + +if __name__ == '__main__': + unittest.main() From 79f7ca7a1e566d8848381a108737276666504105 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Tue, 15 Dec 2015 19:35:18 -0800 Subject: [PATCH 19/30] Adden dpred to be written in the saveDict directive --- SimPEG/Directives.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index e91ffb71..46576df5 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -235,7 +235,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): # Save the file as a npz - np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel) + np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred) From 43424503600e392dd4300a302ffdab65419b9b9b Mon Sep 17 00:00:00 2001 From: GudniRos Date: Thu, 17 Dec 2015 23:02:52 -0800 Subject: [PATCH 20/30] Updated test_MeshIO to remove the temp files after using them. --- tests/mesh/test_MeshIO.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/tests/mesh/test_MeshIO.py b/tests/mesh/test_MeshIO.py index e3e49548..43c0e33e 100644 --- a/tests/mesh/test_MeshIO.py +++ b/tests/mesh/test_MeshIO.py @@ -1,5 +1,5 @@ import numpy as np -import unittest +import unittest, os import SimPEG as simpeg from SimPEG.Mesh import TensorMesh, TreeMesh @@ -29,17 +29,20 @@ class TestOcTreeIO(unittest.TestCase): 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) try: - simpeg.Utils.meshutils.writeVTUFile('test.vtu',mesh,{'arange':vec}) + simpeg.Utils.meshutils.writeVTUFile('temp.vtu',mesh,{'arange':vec}) run = True except: run = False assert run print 'Writing of VTU files is working' + os.remove('temp.vtu') From dedabcc15f7dc97d3f3346045b721669d0335bf9 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 10 Jan 2016 18:36:02 -0800 Subject: [PATCH 21/30] Allow solver kwargs to go to the class directly. --- SimPEG/InvProblem.py | 4 ++-- SimPEG/Utils/SolverUtils.py | 21 +++++++++++++++++---- 2 files changed, 19 insertions(+), 6 deletions(-) diff --git a/SimPEG/InvProblem.py b/SimPEG/InvProblem.py index 32a2195c..0296bf4b 100644 --- a/SimPEG/InvProblem.py +++ b/SimPEG/InvProblem.py @@ -66,8 +66,8 @@ class BaseInvProblem(object): self.curModel = m0 print """SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv. - ***Done using same solver as the problem***""" - self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel)) + ***Done using same Solver and solverOpts as the problem***""" + self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel), **self.prob.solverOpts) @property def warmstart(self): diff --git a/SimPEG/Utils/SolverUtils.py b/SimPEG/Utils/SolverUtils.py index 26ff3e2a..47c899d3 100644 --- a/SimPEG/Utils/SolverUtils.py +++ b/SimPEG/Utils/SolverUtils.py @@ -26,7 +26,14 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6): def __init__(self, A, **kwargs): self.A = A.tocsc() + + self.checkAccuracy = kwargs.get("checkAccuracy", checkAccuracy) + if kwargs.has_key("checkAccuracy"): del kwargs["checkAccuracy"] + self.accuracyTol = kwargs.get("accuracyTol", accuracyTol) + if kwargs.has_key("accuracyTol"): del kwargs["accuracyTol"] + self.kwargs = kwargs + if factorize: self.solver = fun(self.A, **kwargs) @@ -57,8 +64,8 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6): else: X[:,i] = fun(self.A, b[:,i], **self.kwargs) - if checkAccuracy: - _checkAccuracy(self.A, b, X, accuracyTol) + if self.checkAccuracy: + _checkAccuracy(self.A, b, X, self.accuracyTol) return X def clean(self): @@ -81,6 +88,12 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5): def __init__(self, A, **kwargs): self.A = A + + self.checkAccuracy = kwargs.get("checkAccuracy", checkAccuracy) + if kwargs.has_key("checkAccuracy"): del kwargs["checkAccuracy"] + self.accuracyTol = kwargs.get("accuracyTol", accuracyTol) + if kwargs.has_key("accuracyTol"): del kwargs["accuracyTol"] + self.kwargs = kwargs def __mul__(self, b): @@ -108,8 +121,8 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5): else: X[:,i] = out - if checkAccuracy: - _checkAccuracy(self.A, b, X, accuracyTol) + if self.checkAccuracy: + _checkAccuracy(self.A, b, X, self.accuracyTol) return X def clean(self): From 1700f4f9c068464b770f7f9909d6eb20a29331b7 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 14 Jan 2016 14:00:55 -0800 Subject: [PATCH 22/30] add Ainv.clean() to fdem fields, jvec, jtvec --- SimPEG/EM/FDEM/FDEM.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index f2167fd8..3295fd0c 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -50,7 +50,7 @@ class BaseFDEMProblem(BaseEMProblem): Srcs = self.survey.getSrcByFreq(freq) ftype = self._fieldType + 'Solution' F[Srcs, ftype] = sol - + Ainv.clean() return F def Jvec(self, m, v, f=None): @@ -89,6 +89,7 @@ class BaseFDEMProblem(BaseEMProblem): Jv[src, rx] = P(Df_Dm) + Ainv.clean() return Utils.mkvc(Jv) def Jtvec(self, m, v, f=None): @@ -139,7 +140,8 @@ class BaseFDEMProblem(BaseEMProblem): Jtv += - np.array(du_dmT,dtype=complex).real else: raise Exception('Must be real or imag') - + + ATinv.clean() return Jtv def getSourceTerm(self, freq): From 01b1122fcff1be9034e3a5f7ffd7841b6f74859e Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 14 Jan 2016 15:12:09 -0800 Subject: [PATCH 23/30] Add default interpolation location (CC). --- SimPEG/Mesh/TensorMesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index c76306fe..5ea2d86a 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -215,7 +215,7 @@ class BaseTensorMesh(BaseMesh): inside = inside & (pts[:,i] >= tensor.min()-TOL) & (pts[:,i] <= tensor.max()+TOL) return inside - def getInterpolationMat(self, loc, locType, zerosOutside=False): + def getInterpolationMat(self, loc, locType='CC', zerosOutside=False): """ Produces interpolation matrix :param numpy.ndarray loc: Location of points to interpolate to From e15913cf8488daf1dc3c940af903b991ee089bc9 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 14 Jan 2016 15:12:34 -0800 Subject: [PATCH 24/30] Import all code utils into the utils namespace. --- SimPEG/Utils/__init__.py | 2 +- SimPEG/Utils/codeutils.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 250146c1..18c1994f 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, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile +from meshutils import * from curvutils import volTetra, faceInfo, indexCube from interputils import interpmat from CounterUtils import * diff --git a/SimPEG/Utils/codeutils.py b/SimPEG/Utils/codeutils.py index 4a9a28a7..bfd00889 100644 --- a/SimPEG/Utils/codeutils.py +++ b/SimPEG/Utils/codeutils.py @@ -17,7 +17,7 @@ def memProfileWrapper(towrap, *funNames): For example:: - foo_mem = memProfile(foo,'my_func') + foo_mem = memProfileWrapper(foo,['my_func']) fooi = foo_mem() for i in range(5): fooi.my_func() From c6e90230d4b22a615cebfcffee43c624fc646863 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 14 Jan 2016 15:57:15 -0800 Subject: [PATCH 25/30] Updates to the correct pointer for the flow module (docs) --- docs/flow/index.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/flow/index.rst b/docs/flow/index.rst index 3c8083fc..d4299b0d 100644 --- a/docs/flow/index.rst +++ b/docs/flow/index.rst @@ -41,7 +41,7 @@ Here we reproduce the results from Celia et al. (1990): Richards ======== -.. automodule:: simpegFLOW.Richards.Empirical +.. automodule:: SimPEG.FLOW.Richards.Empirical :show-inheritance: :members: :undoc-members: From 09161ff68e657e95a66dc283914f309ad903eac7 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 28 Jan 2016 13:58:32 -0800 Subject: [PATCH 26/30] Hopefully get a better error message on travis. --- tests/mesh/test_MeshIO.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/tests/mesh/test_MeshIO.py b/tests/mesh/test_MeshIO.py index 43c0e33e..52e5740d 100644 --- a/tests/mesh/test_MeshIO.py +++ b/tests/mesh/test_MeshIO.py @@ -35,12 +35,7 @@ class TestOcTreeIO(unittest.TestCase): def test_VTUfiles(self): mesh = self.mesh vec = np.arange(mesh.nC) - try: - simpeg.Utils.meshutils.writeVTUFile('temp.vtu',mesh,{'arange':vec}) - run = True - except: - run = False - assert run + simpeg.Utils.meshutils.writeVTUFile('temp.vtu',mesh,{'arange':vec}) print 'Writing of VTU files is working' os.remove('temp.vtu') From 860bd5638ab056f90c3b2a85c9165f4625e667af Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 28 Jan 2016 14:19:29 -0800 Subject: [PATCH 27/30] Add VTK to the travis dependencies. I am not adding this to the requirements.txt file. --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index c2c672bf..de145ecb 100644 --- a/.travis.yml +++ b/.travis.yml @@ -26,7 +26,7 @@ before_install: # Install packages install: - - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython + - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython vtk - pip install nose-cov python-coveralls # - pip install -r requirements.txt - python setup.py install From 1bcb572c459c215b451bbff9127af26d7f8be9f0 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 28 Jan 2016 14:29:45 -0800 Subject: [PATCH 28/30] Remove ifmain from TreeMesh --- .gitignore | 1 - SimPEG/Mesh/TreeMesh.py | 125 ---------------------------------------- 2 files changed, 126 deletions(-) diff --git a/.gitignore b/.gitignore index 6a826230..fc798cea 100644 --- a/.gitignore +++ b/.gitignore @@ -39,4 +39,3 @@ nosetests.xml *.sublime-workspace docs/_build/ Makefile -diff.temp diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index a4c4e845..0205e972 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -2335,128 +2335,3 @@ class NotBalancedException(TreeException): pass class CellLookUpException(TreeException): pass - -if __name__ == '__main__': - - - import matplotlib.pyplot as plt - import matplotlib - from mpl_toolkits.mplot3d import Axes3D - import matplotlib.colors as colors - import matplotlib.cm as cmx - - def topo(x): - return np.sin(x*(2.*np.pi))*0.3 + 0.5 - - def function(cell): - r = cell.center - np.array([0.5]*len(cell.center)) - dist = np.sqrt(r.dot(r)) - # dist2 = np.abs(cell.center[-1] - topo(cell.center[0])) - - # dist = min([dist1,dist2]) - # if dist < 0.05: - # return 5 - if dist < 0.1: - return 5 - if dist < 0.2: - return 4 - if dist < 0.4: - return 3 - return 2 - - # T = TreeMesh([[(1,128)],[(1,128)],[(1,128)]],levels=7) - # T = TreeMesh([128,128,128]) - # T = TreeMesh([64,64],levels=6) - T = TreeMesh([8,8]) - # T = TreeMesh([[(1,128)],[(1,128)]],levels=7) - # T.refine(lambda xc:2, balance=False) - # T._index([0,0,0]) - # T._pointer(0) - - - # tic = time.time() - T.refine(function)#, balance=False) - # print time.time() - tic - # print T.nC - # T.plotSlice(np.log(T.vol))#np.random.rand(T.nC)) - T.plotGrid() - # print [c for c in T] - c = T[0] - plt.plot(c.center[0],c.center[1],'r.') - nodes = c.nodes - for n in nodes: - _ = T._gridN[n,:] - plt.plot(_[0],_[1],'gs') - plt.show() - blah - - # T.plotImage(np.arange(len(T.vol)),showIt=True) - - # print T.getFaceInnerProduct() - # print T.gridFz - - - # T._refineCell([8,0,1]) - # T._refineCell([8,0,2]) - # T._refineCell([12,0,2]) - # T._refineCell([8,4,2]) - # T._refineCell([6,0,3]) - # T._refineCell([8,8,1]) - # T._refineCell([0,0,0,1]) - # T.__dirty__ = True - - - # print T.gridFx.shape[0], T.nFx - - - - ax = plt.subplot(211) - ax.spy(T.edgeCurl) - - # print Mesh.TensorMesh([2,2,2]).edgeCurl.todense() - # print T.edgeCurl.todense() - # print Mesh.TensorMesh([2,2,2]).edgeCurl.todense() - T.edgeCurl.todense() - # print T.gridEy - Mesh.TensorMesh([2,2,2]).gridEy - - # print T.edge - # T.plotGrid(ax=ax) - - # R = deflationMatrix(T._facesX, T._hangingFx, T._fx2i) - # print R - - ax = plt.subplot(212)#, projection='3d') - ax.spy(Mesh.TensorMesh([2,2,2]).edgeCurl) - - # ax = plt.subplot(313) - # ax.spy(T.faceDiv[:,:T.nFx] * R) - - - # T.balance() - # T.plotGrid(ax=ax) - - # cx = T._getNextCell([0,0,1],direction=0,positive=True) - # print cx - # # print [T._asPointer(_) for _ in cx] - # cx = T._getNextCell([8,0,3],direction=0,positive=False) - # print T._asPointer(cx) - # cx = T._getNextCell([8,8,1],direction=1,positive=False) - # print cx, #[T._asPointer(_) for _ in cx] - # cm = T._getNextCell([64,80,4],direction=0,positive=False) - # cy = T._getNextCell([64,80,4],direction=1,positive=True) - # cp = T._getNextCell([64,80,4],direction=1,positive=False) - - # ax.plot( T._cellN([4,0,1])[0],T._cellN([4,0,1])[1], 'yd') - # ax.plot( T._cellN(cx)[0],T._cellN(cx)[1], 'ys') - # ax.plot( T._cellN(cm)[0],T._cellN(cm)[1], 'ys') - # ax.plot( T._cellN(cy)[0],T._cellN(cy)[1], 'ys') - # ax.plot( T._cellN(cp[0])[0],T._cellN(cp[0])[1], 'ys') - # ax.plot( T._cellN(cp[1])[0],T._cellN(cp[1])[1], 'ys') - - - - - - # print T.nN - - plt.show() - From 43c49d5f157dc7ea803fd6a420a0c38cb4449c00 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 28 Jan 2016 17:53:10 -0800 Subject: [PATCH 29/30] Address mesh IO #212 --- SimPEG/Mesh/MeshIO.py | 416 ++++++++++++++ SimPEG/Mesh/TensorMesh.py | 1117 +++++++++++++++++++------------------ SimPEG/Mesh/TreeMesh.py | 3 +- SimPEG/Utils/meshutils.py | 402 ------------- tests/mesh/test_MeshIO.py | 67 ++- 5 files changed, 1038 insertions(+), 967 deletions(-) create mode 100644 SimPEG/Mesh/MeshIO.py diff --git a/SimPEG/Mesh/MeshIO.py b/SimPEG/Mesh/MeshIO.py new file mode 100644 index 00000000..7501a66f --- /dev/null +++ b/SimPEG/Mesh/MeshIO.py @@ -0,0 +1,416 @@ +import numpy as np, os +from SimPEG import Utils + +class TensorMeshIO(object): + + @classmethod + def readUBC(TensorMesh, fileName): + """ + Read UBC GIF 3DTensor mesh and generate 3D Tensor mesh in simpegTD + + Input: + :param fileName, path to the UBC GIF mesh file + + Output: + :param SimPEG TensorMesh object + """ + + # Interal function to read cell size lines for the UBC mesh files. + def readCellLine(line): + for seg in line.split(): + if '*' in seg: + st = seg + sp = seg.split('*') + re = np.array(sp[0],dtype=int)*(' ' + sp[1]) + line = line.replace(st,re.strip()) + return np.array(line.split(),dtype=float) + + # Read the file as line strings, remove lines with comment = ! + msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!') + + # Fist line is the size of the model + sizeM = np.array(msh[0].split(),dtype=float) + # Second line is the South-West-Top corner coordinates. + x0 = np.array(msh[1].split(),dtype=float) + # Read the cell sizes + h1 = readCellLine(msh[2]) + h2 = readCellLine(msh[3]) + h3temp = readCellLine(msh[4]) + h3 = h3temp[::-1] # Invert the indexing of the vector to start from the bottom. + # Adjust the reference point to the bottom south west corner + x0[2] = x0[2] - np.sum(h3) + # Make the mesh + tensMsh = TensorMesh([h1,h2,h3],x0) + return tensMsh + + @classmethod + def readVTK(TensorMesh, fileName): + """ + Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model + + Input: + :param vtrFileName, path to the vtr model file to write to + + Output: + :return SimPEG TensorMesh object + :return SimPEG model dictionary + + """ + # Import + from vtk import vtkXMLRectilinearGridReader as vtrFileReader + from vtk.util.numpy_support import vtk_to_numpy + + # Read the file + vtrReader = vtrFileReader() + vtrReader.SetFileName(fileName) + vtrReader.Update() + vtrGrid = vtrReader.GetOutput() + # Sort information + hx = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetXCoordinates()))) + xR = vtk_to_numpy(vtrGrid.GetXCoordinates())[0] + hy = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetYCoordinates()))) + yR = vtk_to_numpy(vtrGrid.GetYCoordinates())[0] + zD = np.diff(vtk_to_numpy(vtrGrid.GetZCoordinates())) + # Check the direction of hz + if np.all(zD < 0): + hz = np.abs(zD[::-1]) + zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[-1] + else: + hz = np.abs(zD) + zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[0] + x0 = np.array([xR,yR,zR]) + + # Make the SimPEG object + tensMsh = TensorMesh([hx,hy,hz],x0) + + # Grap the models + models = {} + for i in np.arange(vtrGrid.GetCellData().GetNumberOfArrays()): + modelName = vtrGrid.GetCellData().GetArrayName(i) + if np.all(zD < 0): + modFlip = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i)) + tM = tensMsh.r(modFlip,'CC','CC','M') + modArr = tensMsh.r(tM[:,:,::-1],'CC','CC','V') + else: + modArr = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i)) + models[modelName] = modArr + + # Return the data + return tensMsh, models + + def writeVTK(mesh, fileName, models=None): + """ + Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model. + + Input: + :param str, path to the output vtk file + :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK + :param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells + + """ + # Import + from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter, VTK_VERSION + from vtk.util.numpy_support import numpy_to_vtk + + # Deal with dimensionalities + if mesh.dim >= 1: + vX = mesh.vectorNx + xD = mesh.nNx + yD,zD = 1,1 + vY, vZ = np.array([0,0]) + if mesh.dim >= 2: + vY = mesh.vectorNy + yD = mesh.nNy + if mesh.dim == 3: + vZ = mesh.vectorNz + zD = mesh.nNz + # Use rectilinear VTK grid. + # Assign the spatial information. + vtkObj = rectGrid() + vtkObj.SetDimensions(xD,yD,zD) + vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1)) + vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1)) + vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1)) + + # Assign the model('s) to the object + if models is not None: + for item in models.iteritems(): + # Convert numpy array + vtkDoubleArr = numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtkObj.GetCellData().AddArray(vtkDoubleArr) + # Set the active scalar + vtkObj.GetCellData().SetActiveScalars(models.keys()[0]) + # vtkObj.Update() + + # Check the extension of the fileName + ext = os.path.splitext(fileName)[1] + if ext is '': + fileName = fileName + '.vtr' + elif ext not in '.vtr': + raise IOError('{:s} is an incorrect extension, has to be .vtr') + # Write the file. + vtrWriteFilter = rectWriter() + if float(VTK_VERSION.split('.')[0]) >=6: + vtrWriteFilter.SetInputData(vtkObj) + else: + vtuWriteFilter.SetInput(vtuObj) + vtrWriteFilter.SetFileName(fileName) + vtrWriteFilter.Update() + + + def readModelUBC(mesh, fileName): + """ + Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg + + Input: + :param fileName, path to the UBC GIF mesh file to read + :param mesh, TensorMesh object, mesh that coresponds to the model + + Output: + :return numpy array, model with TensorMesh ordered + """ + f = open(fileName, 'r') + model = np.array(map(float, f.readlines())) + f.close() + model = np.reshape(model, (mesh.nCz, mesh.nCx, mesh.nCy), order = 'F') + model = model[::-1,:,:] + model = np.transpose(model, (1, 2, 0)) + model = Utils.mkvc(model) + return model + + def writeModelUBC(mesh, fileName, model): + """ + Writes a model associated with a SimPEG TensorMesh + to a UBC-GIF format model file. + + :param str fileName: File to write to + :param simpeg.Mesh.TensorMesh mesh: The mesh + :param numpy.ndarray model: The model + """ + + # Reshape model to a matrix + modelMat = mesh.r(model,'CC','CC','M') + # Transpose the axes + modelMatT = modelMat.transpose((2,0,1)) + # Flip z to positive down + modelMatTR = Utils.mkvc(modelMatT[::-1,:,:]) + + np.savetxt(fileName, modelMatTR.ravel()) + + def writeUBC(mesh, fileName, models=None): + """ + Writes a SimPEG TensorMesh to a UBC-GIF format mesh file. + + :param str fileName: File to write to + :param simpeg.Mesh.TensorMesh mesh: The mesh + + """ + assert mesh.dim == 3 + s = '' + s += '%i %i %i\n' %tuple(mesh.vnC) + origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated. + origin.dtype = float + + 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[::-1]) + f = open(fileName, 'w') + f.write(s) + f.close() + + if models is None: return + assert type(models) is dict, 'models must be a dict' + for key in models: + assert type(key) is str, 'The dict key is a file name' + mesh.writeModelUBC(key, models[key]) + +class TreeMeshIO(object): + + def writeUBC(mesh, fileName, models=None): + """ + Write UBC ocTree mesh and model files from a simpeg ocTree mesh and model. + + :param str fileName: File to write to + :param simpeg.Mesh.TreeMesh mesh: The mesh + :param dictionary models: The models in a dictionary, where the keys is the name of the of the model file + """ + + # Calculate information to write in the file. + # Number of cells in the underlying mesh + nCunderMesh = np.array([h.size for h in mesh.h],dtype=np.int64) + # The top-south-west most corner of the mesh + tswCorn = mesh.x0 + np.array([0,0,np.sum(mesh.h[2])]) + # Smallest cell size + smallCell = np.array([h.min() for h in mesh.h]) + # Number of cells + nrCells = mesh.nC + + ## Extract iformation about the cells. + # cell pointers + cellPointers = np.array([c._pointer for c in mesh]) + # cell with + cellW = np.array([ mesh._levelWidth(i) for i in cellPointers[:,-1] ]) + # Need to shift the pointers to work with UBC indexing + # UBC Octree indexes always the top-left-close (top-south-west) corner first and orders the cells in z(top-down),x,y vs x,y,z(bottom-up). + # Shift index up by 1 + ubcCellPt = cellPointers[:,0:-1].copy() + np.array([1.,1.,1.]) + # Need reindex the z index to be from the top-left-close corner and to be from the global top. + ubcCellPt[:,2] = ( nCunderMesh[-1] + 2) - (ubcCellPt[:,2] + cellW) + + # Reorder the ubcCellPt + ubcReorder = np.argsort(ubcCellPt.view(','.join(3*['float'])),axis=0,order=['f2','f1','f0'])[:,0] + # Make a array with the pointers and the withs, that are order in the ubc ordering + indArr = np.concatenate((ubcCellPt[ubcReorder,:],cellW[ubcReorder].reshape((-1,1)) ),axis=1) + + ## Write the UBC octree mesh file + with open(fileName,'w') as mshOut: + mshOut.write('{:.0f} {:.0f} {:.0f}\n'.format(nCunderMesh[0],nCunderMesh[1],nCunderMesh[2])) + mshOut.write('{:.4f} {:.4f} {:.4f}\n'.format(tswCorn[0],tswCorn[1],tswCorn[2])) + mshOut.write('{:.3f} {:.3f} {:.3f}\n'.format(smallCell[0],smallCell[1],smallCell[2])) + mshOut.write('{:.0f} \n'.format(nrCells)) + np.savetxt(mshOut,indArr,fmt='%i') + + ## Print the models + # Assign the model('s) to the object + if models is not None: + # indUBCvector = np.argsort(cX0[np.argsort(np.concatenate((cX0[:,0:2],cX0[:,2:3].max() - cX0[:,2:3]),axis=1).view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0]].view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0] + for item in models.iteritems(): + # Save the data + np.savetxt(item[0],item[1][ubcReorder],fmt='%3.5e') + + @classmethod + def readUBC(TreeMesh, meshFile): + """ + Read UBC 3D OcTree mesh and/or modelFiles + + Input: + :param str meshFile: path to the UBC GIF OcTree mesh file to read + + Output: + :return SimPEG.Mesh.TreeMesh mesh: The octree mesh + :return list of ndarray's: models as a list of numpy array's + """ + + ## Read the file lines + fileLines = np.genfromtxt(meshFile,dtype=str,delimiter='\n') + # Extract the data + nCunderMesh = np.array(fileLines[0].split(),dtype=float) + # I think this is the case? + if np.unique(nCunderMesh).size >1: + raise Exception('SimPEG TreeMeshes have the same number of cell in all directions') + tswCorn = np.array(fileLines[1].split(),dtype=float) + smallCell = np.array(fileLines[2].split(),dtype=float) + nrCells = np.array(fileLines[3].split(),dtype=float) + # Read the index array + indArr = np.genfromtxt(fileLines[4::],dtype=np.int) + + ## Calculate simpeg parameters + h1,h2,h3 = [np.ones(nr)*sz for nr,sz in zip(nCunderMesh,smallCell)] + x0 = tswCorn - np.array([0,0,np.sum(h3)]) + # Need to convert the index array to a points list that complies with SimPEG TreeMesh. + # Shift to start at 0 + simpegCellPt = indArr[:,0:-1].copy() + simpegCellPt[:,2] = ( nCunderMesh[-1] + 2) - (simpegCellPt[:,2] + indArr[:,3]) + # Need reindex the z index to be from the bottom-left-close corner and to be from the global bottom. + simpegCellPt = simpegCellPt - np.array([1.,1.,1.]) + + # Calculate the cell level + simpegLevel = np.log2(np.min(nCunderMesh)) - np.log2(indArr[:,3]) + # Make a pointer matrix + simpegPointers = np.concatenate((simpegCellPt,simpegLevel.reshape((-1,1))),axis=1) + + ## Make the tree mesh + mesh = TreeMesh([h1,h2,h3],x0) + mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()]) + + # Figure out the reordering + mesh._simpegReorderUBC = np.argsort(np.array([mesh._index(i) for i in simpegPointers.tolist()])) + # mesh._simpegReorderUBC = np.argsort((np.array([[1,1,1,-1]])*simpegPointers).view(','.join(4*['float'])),axis=0,order=['f3','f2','f1','f0'])[:,0] + + return mesh + + + def readModelUBC(mesh, fileName): + """ + Read UBC OcTree model and get vector + + Input: + :param fileName, path to the UBC GIF model file to read + + Output: + :return numpy array, OcTree model + """ + + if type(fileName) is list: + out = {} + for f in fileName: + out[f] = mesh.readModelUBC(f) + return out + + assert hasattr(mesh, '_simpegReorderUBC'), 'The file must have been loaded from a UBC format.' + assert mesh.dim == 3 + + modList = [] + modArr = np.loadtxt(fileName) + if len(modArr.shape) == 1: + modList.append(modArr[mesh._simpegReorderUBC]) + else: + modList.append(modArr[mesh._simpegReorderUBC,:]) + return modList + + def writeVTK(mesh, fileName, models=None): + """ + Function to write a VTU file from a SimPEG TreeMesh and model. + """ + import vtk + from vtk import vtkXMLUnstructuredGridWriter as Writer, VTK_VERSION + from vtk.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray + + if str(type(mesh)).split()[-1][1:-2] not in 'SimPEG.Mesh.TreeMesh.TreeMesh': + raise IOError('mesh is not a SimPEG TreeMesh.') + + # Make the data parts for the vtu object + # Points + mesh.number() + ptsMat = mesh._gridN + mesh.x0 + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(numpy_to_vtk(ptsMat,deep=True)) + # Cells + cellConn = np.array([c.nodes for c in mesh],dtype=np.int64) + + cellsMat = np.concatenate((np.ones((cellConn.shape[0],1),dtype=np.int64)*cellConn.shape[1],cellConn),axis=1).ravel() + cellsArr = vtk.vtkCellArray() + cellsArr.SetNumberOfCells(cellConn.shape[0]) + cellsArr.SetCells(cellConn.shape[0],numpy_to_vtkIdTypeArray(cellsMat,deep=True)) + + # Make the object + vtuObj = vtk.vtkUnstructuredGrid() + vtuObj.SetPoints(vtkPts) + vtuObj.SetCells(vtk.VTK_VOXEL,cellsArr) + # Add the level of refinement as a cell array + cellSides = np.array([np.array(vtuObj.GetCell(i).GetBounds()).reshape((3,2)).dot(np.array([-1, 1])) for i in np.arange(vtuObj.GetNumberOfCells())]) + uniqueLevel, indLevel = np.unique(np.prod(cellSides,axis=1),return_inverse=True) + refineLevelArr = numpy_to_vtk(indLevel.max() - indLevel,deep=1) + refineLevelArr.SetName('octreeLevel') + vtuObj.GetCellData().AddArray(refineLevelArr) + # Assign the model('s) to the object + if models is not None: + for item in models.iteritems(): + # Convert numpy array + vtkDoubleArr = numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtuObj.GetCellData().AddArray(vtkDoubleArr) + + # Make the writer + vtuWriteFilter = Writer() + if float(VTK_VERSION.split('.')[0]) >=6: + vtuWriteFilter.SetInputData(vtuObj) + else: + vtuWriteFilter.SetInput(vtuObj) + vtuWriteFilter.SetFileName(fileName) + # Write the file + vtuWriteFilter.Update() + diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 5ea2d86a..508f015c 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -1,558 +1,559 @@ -from SimPEG import Utils, np, sp -from BaseMesh import BaseMesh, BaseRectangularMesh -from View import TensorView -from DiffOperators import DiffOperators -from InnerProducts import InnerProducts - -class BaseTensorMesh(BaseMesh): - - __metaclass__ = Utils.SimPEGMetaClass - - _meshType = 'BASETENSOR' - - _unitDimensions = [1, 1, 1] - - def __init__(self, h_in, x0_in=None): - assert type(h_in) in [list, tuple], '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 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 isinstance(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) - - if isinstance(self, BaseRectangularMesh): - BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0) - else: - BaseMesh.__init__(self, np.array([x.size for x in h]), x0) - - # Ensure h contains 1D vectors - self._h = [Utils.mkvc(x.astype(float)) for x in h] - - @property - def h(self): - """h is a list containing the cell widths of the tensor mesh in each dimension.""" - return self._h - - @property - def hx(self): - "Width of cells in the x direction" - return self._h[0] - - @property - def hy(self): - "Width of cells in the y direction" - return None if self.dim < 2 else self._h[1] - - @property - def hz(self): - "Width of cells in the z direction" - return None if self.dim < 3 else self._h[2] - - @property - def vectorNx(self): - """Nodal grid vector (1D) in the x direction.""" - return np.r_[0., self.hx.cumsum()] + self.x0[0] - - @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] - - @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] - - @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] - - @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] - - @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] - - @property - def gridCC(self): - """Cell-centered grid.""" - return self._getTensorGrid('CC') - - @property - def gridN(self): - """Nodal grid.""" - return self._getTensorGrid('N') - - @property - def gridFx(self): - """Face staggered grid in the x direction.""" - if self.nFx == 0: return - return self._getTensorGrid('Fx') - - @property - def gridFy(self): - """Face staggered grid in the y direction.""" - if self.nFy == 0 or self.dim < 2: return - return self._getTensorGrid('Fy') - - @property - def gridFz(self): - """Face staggered grid in the z direction.""" - if self.nFz == 0 or self.dim < 3: return - return self._getTensorGrid('Fz') - - @property - def gridEx(self): - """Edge staggered grid in the x direction.""" - if self.nEx == 0: return - return self._getTensorGrid('Ex') - - @property - def gridEy(self): - """Edge staggered grid in the y direction.""" - if self.nEy == 0 or self.dim < 2: return - return self._getTensorGrid('Ey') - - @property - def gridEz(self): - """Edge staggered grid in the z direction.""" - if self.nEz == 0 or self.dim < 3: return - return self._getTensorGrid('Ez') - - 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) - - def getTensor(self, key): - """ Returns a tensor list. - - :param str key: What tensor (see below) - :rtype: list - :return: list of the tensors that make up the mesh. - - key can be:: - - 'CC' -> scalar field defined on cell centers - 'N' -> scalar field defined on nodes - 'Fx' -> x-component of field defined on faces - 'Fy' -> y-component of field defined on faces - 'Fz' -> z-component of field defined on faces - 'Ex' -> x-component of field defined on edges - 'Ey' -> y-component of field defined on edges - 'Ez' -> z-component of field defined on edges - - """ - - if key == 'Fx': - ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] - elif key == 'Fy': - ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] - elif key == 'Fz': - ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] - elif key == 'Ex': - ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] - elif key == 'Ey': - ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] - elif key == 'Ez': - ten = [self.vectorNx , self.vectorNy , self.vectorCCz] - elif key == 'CC': - ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] - elif key == 'N': - ten = [self.vectorNx , self.vectorNy , self.vectorNz ] - - 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 - """ - pts = Utils.asArray_N_x_Dim(pts, self.dim) - - tensors = self.getTensor(locType) - - if locType == 'N' and self._meshType == 'CYL': - #NOTE: for a CYL mesh we add a node to check if we are inside in the radial direction! - tensors[0] = np.r_[0.,tensors[0]] - tensors[1] = np.r_[tensors[1], 2.0*np.pi] - - inside = np.ones(pts.shape[0],dtype=bool) - for i, tensor in enumerate(tensors): - TOL = np.diff(tensor).min() * 1.0e-10 - inside = inside & (pts[:,i] >= tensor.min()-TOL) & (pts[:,i] <= tensor.max()+TOL) - return inside - - def getInterpolationMat(self, loc, locType='CC', 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 self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']: - raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType) - - loc = Utils.asArray_N_x_Dim(loc, self.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)) - # 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() - - - def _fastInnerProduct(self, projType, prop=None, invProp=False, invMat=False): - """ - Fast version of getFaceInnerProduct. - This does not handle the case of a full tensor prop. - - :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param str projType: 'E' or 'F' - :param bool returnP: returns the projection matrices - :param bool invProp: inverts the material property - :param bool invMat: inverts the matrix - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ - assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" - - if prop is None: - prop = np.ones(self.nC) - - if invProp: - prop = 1./prop - - if Utils.isScalar(prop): - prop = prop*np.ones(self.nC) - - if prop.size == self.nC: - Av = getattr(self, 'ave'+projType+'2CC') - Vprop = self.vol * Utils.mkvc(prop) - M = self.dim * Utils.sdiag(Av.T * Vprop) - elif prop.size == self.nC*self.dim: - Av = getattr(self, 'ave'+projType+'2CCV') - V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - M = Utils.sdiag(Av.T * V * Utils.mkvc(prop)) - else: - return None - - if invMat: - return Utils.sdInv(M) - else: - return M - - def _fastInnerProductDeriv(self, projType, prop, invProp=False, invMat=False): - """ - :param str projType: 'E' or 'F' - :param TensorType tensorType: type of the tensor - :param bool invProp: inverts the material property - :param bool invMat: inverts the matrix - :rtype: function - :return: dMdmu, the derivative of the inner product matrix - """ - assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" - tensorType = Utils.TensorType(self, prop) - - dMdprop = None - - if invMat: - MI = self._fastInnerProduct(projType, prop, invProp=invProp, invMat=invMat) - - if tensorType == 0: - Av = getattr(self, 'ave'+projType+'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 not invMat and not invProp: - dMdprop = self.dim * Av.T * V * ones - elif invMat and invProp: - dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * ones * Utils.sdiag(1./prop**2) - - if tensorType == 1: - Av = getattr(self, 'ave'+projType+'2CC') - V = Utils.sdiag(self.vol) - if not invMat and not invProp: - dMdprop = self.dim * Av.T * V - elif invMat and invProp: - dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2) - - if tensorType == 2: # anisotropic - Av = getattr(self, 'ave'+projType+'2CCV') - V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - if not invMat and not invProp: - dMdprop = Av.T * V - elif invMat and invProp: - dMdprop = Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2) - - if dMdprop is not None: - def innerProductDeriv(v=None): - if v is None: - print 'Depreciation Warning: TensorMesh.innerProductDeriv. You should be supplying a vector. Use: sdiag(u)*dMdprop' - return dMdprop - return Utils.sdiag(v) * dMdprop - return innerProductDeriv - else: - return None - - - -class TensorMesh(BaseTensorMesh, 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 using :func:`SimPEG.Utils.meshutils.meshTensor`: - - .. plot:: - :include-source: - - from SimPEG import Mesh, Utils - 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:: - - 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 += ' {0:.2f},'.format(h) - else: - outStr += ' {0:d}*{1:.2f},'.format(n,h) - - return outStr[:-1] - - if self.dim == 1: - outStr += '\n x0: {0:.2f}'.format(self.x0[0]) - outStr += '\n nCx: {0:d}'.format(self.nCx) - outStr += printH(self.hx, outStr='\n hx:') - pass - elif self.dim == 2: - outStr += '\n x0: {0:.2f}'.format(self.x0[0]) - outStr += '\n y0: {0:.2f}'.format(self.x0[1]) - outStr += '\n nCx: {0:d}'.format(self.nCx) - outStr += '\n nCy: {0:d}'.format(self.nCy) - outStr += printH(self.hx, outStr='\n hx:') - outStr += printH(self.hy, outStr='\n hy:') - elif self.dim == 3: - outStr += '\n x0: {0:.2f}'.format(self.x0[0]) - outStr += '\n y0: {0:.2f}'.format(self.x0[1]) - outStr += '\n z0: {0:.2f}'.format(self.x0[2]) - outStr += '\n nCx: {0:d}'.format(self.nCx) - outStr += '\n nCy: {0:d}'.format(self.nCy) - outStr += '\n nCz: {0:d}'.format(self.nCz) - outStr += printH(self.hx, outStr='\n hx:') - outStr += printH(self.hy, outStr='\n hy:') - 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 +from SimPEG import Utils, np, sp +from BaseMesh import BaseMesh, BaseRectangularMesh +from View import TensorView +from DiffOperators import DiffOperators +from InnerProducts import InnerProducts +from MeshIO import TensorMeshIO + +class BaseTensorMesh(BaseMesh): + + __metaclass__ = Utils.SimPEGMetaClass + + _meshType = 'BASETENSOR' + + _unitDimensions = [1, 1, 1] + + def __init__(self, h_in, x0_in=None): + assert type(h_in) in [list, tuple], '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 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 isinstance(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) + + if isinstance(self, BaseRectangularMesh): + BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0) + else: + BaseMesh.__init__(self, np.array([x.size for x in h]), x0) + + # Ensure h contains 1D vectors + self._h = [Utils.mkvc(x.astype(float)) for x in h] + + @property + def h(self): + """h is a list containing the cell widths of the tensor mesh in each dimension.""" + return self._h + + @property + def hx(self): + "Width of cells in the x direction" + return self._h[0] + + @property + def hy(self): + "Width of cells in the y direction" + return None if self.dim < 2 else self._h[1] + + @property + def hz(self): + "Width of cells in the z direction" + return None if self.dim < 3 else self._h[2] + + @property + def vectorNx(self): + """Nodal grid vector (1D) in the x direction.""" + return np.r_[0., self.hx.cumsum()] + self.x0[0] + + @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] + + @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] + + @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] + + @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] + + @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] + + @property + def gridCC(self): + """Cell-centered grid.""" + return self._getTensorGrid('CC') + + @property + def gridN(self): + """Nodal grid.""" + return self._getTensorGrid('N') + + @property + def gridFx(self): + """Face staggered grid in the x direction.""" + if self.nFx == 0: return + return self._getTensorGrid('Fx') + + @property + def gridFy(self): + """Face staggered grid in the y direction.""" + if self.nFy == 0 or self.dim < 2: return + return self._getTensorGrid('Fy') + + @property + def gridFz(self): + """Face staggered grid in the z direction.""" + if self.nFz == 0 or self.dim < 3: return + return self._getTensorGrid('Fz') + + @property + def gridEx(self): + """Edge staggered grid in the x direction.""" + if self.nEx == 0: return + return self._getTensorGrid('Ex') + + @property + def gridEy(self): + """Edge staggered grid in the y direction.""" + if self.nEy == 0 or self.dim < 2: return + return self._getTensorGrid('Ey') + + @property + def gridEz(self): + """Edge staggered grid in the z direction.""" + if self.nEz == 0 or self.dim < 3: return + return self._getTensorGrid('Ez') + + 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) + + def getTensor(self, key): + """ Returns a tensor list. + + :param str key: What tensor (see below) + :rtype: list + :return: list of the tensors that make up the mesh. + + key can be:: + + 'CC' -> scalar field defined on cell centers + 'N' -> scalar field defined on nodes + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces + 'Ex' -> x-component of field defined on edges + 'Ey' -> y-component of field defined on edges + 'Ez' -> z-component of field defined on edges + + """ + + if key == 'Fx': + ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] + elif key == 'Fy': + ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] + elif key == 'Fz': + ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] + elif key == 'Ex': + ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] + elif key == 'Ey': + ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] + elif key == 'Ez': + ten = [self.vectorNx , self.vectorNy , self.vectorCCz] + elif key == 'CC': + ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] + elif key == 'N': + ten = [self.vectorNx , self.vectorNy , self.vectorNz ] + + 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 + """ + pts = Utils.asArray_N_x_Dim(pts, self.dim) + + tensors = self.getTensor(locType) + + if locType == 'N' and self._meshType == 'CYL': + #NOTE: for a CYL mesh we add a node to check if we are inside in the radial direction! + tensors[0] = np.r_[0.,tensors[0]] + tensors[1] = np.r_[tensors[1], 2.0*np.pi] + + inside = np.ones(pts.shape[0],dtype=bool) + for i, tensor in enumerate(tensors): + TOL = np.diff(tensor).min() * 1.0e-10 + inside = inside & (pts[:,i] >= tensor.min()-TOL) & (pts[:,i] <= tensor.max()+TOL) + return inside + + def getInterpolationMat(self, loc, locType='CC', 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 self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']: + raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType) + + loc = Utils.asArray_N_x_Dim(loc, self.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)) + # 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() + + + def _fastInnerProduct(self, projType, prop=None, invProp=False, invMat=False): + """ + Fast version of getFaceInnerProduct. + This does not handle the case of a full tensor prop. + + :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param str projType: 'E' or 'F' + :param bool returnP: returns the projection matrices + :param bool invProp: inverts the material property + :param bool invMat: inverts the matrix + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" + + if prop is None: + prop = np.ones(self.nC) + + if invProp: + prop = 1./prop + + if Utils.isScalar(prop): + prop = prop*np.ones(self.nC) + + if prop.size == self.nC: + Av = getattr(self, 'ave'+projType+'2CC') + Vprop = self.vol * Utils.mkvc(prop) + M = self.dim * Utils.sdiag(Av.T * Vprop) + elif prop.size == self.nC*self.dim: + Av = getattr(self, 'ave'+projType+'2CCV') + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + M = Utils.sdiag(Av.T * V * Utils.mkvc(prop)) + else: + return None + + if invMat: + return Utils.sdInv(M) + else: + return M + + def _fastInnerProductDeriv(self, projType, prop, invProp=False, invMat=False): + """ + :param str projType: 'E' or 'F' + :param TensorType tensorType: type of the tensor + :param bool invProp: inverts the material property + :param bool invMat: inverts the matrix + :rtype: function + :return: dMdmu, the derivative of the inner product matrix + """ + assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" + tensorType = Utils.TensorType(self, prop) + + dMdprop = None + + if invMat: + MI = self._fastInnerProduct(projType, prop, invProp=invProp, invMat=invMat) + + if tensorType == 0: + Av = getattr(self, 'ave'+projType+'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 not invMat and not invProp: + dMdprop = self.dim * Av.T * V * ones + elif invMat and invProp: + dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * ones * Utils.sdiag(1./prop**2) + + if tensorType == 1: + Av = getattr(self, 'ave'+projType+'2CC') + V = Utils.sdiag(self.vol) + if not invMat and not invProp: + dMdprop = self.dim * Av.T * V + elif invMat and invProp: + dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2) + + if tensorType == 2: # anisotropic + Av = getattr(self, 'ave'+projType+'2CCV') + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + if not invMat and not invProp: + dMdprop = Av.T * V + elif invMat and invProp: + dMdprop = Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2) + + if dMdprop is not None: + def innerProductDeriv(v=None): + if v is None: + print 'Depreciation Warning: TensorMesh.innerProductDeriv. You should be supplying a vector. Use: sdiag(u)*dMdprop' + return dMdprop + return Utils.sdiag(v) * dMdprop + return innerProductDeriv + else: + return None + + + +class TensorMesh(BaseTensorMesh, BaseRectangularMesh, TensorView, DiffOperators, InnerProducts, TensorMeshIO): + """ + 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 using :func:`SimPEG.Utils.meshutils.meshTensor`: + + .. plot:: + :include-source: + + from SimPEG import Mesh, Utils + 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:: + + 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 += ' {0:.2f},'.format(h) + else: + outStr += ' {0:d}*{1:.2f},'.format(n,h) + + return outStr[:-1] + + if self.dim == 1: + outStr += '\n x0: {0:.2f}'.format(self.x0[0]) + outStr += '\n nCx: {0:d}'.format(self.nCx) + outStr += printH(self.hx, outStr='\n hx:') + pass + elif self.dim == 2: + outStr += '\n x0: {0:.2f}'.format(self.x0[0]) + outStr += '\n y0: {0:.2f}'.format(self.x0[1]) + outStr += '\n nCx: {0:d}'.format(self.nCx) + outStr += '\n nCy: {0:d}'.format(self.nCy) + outStr += printH(self.hx, outStr='\n hx:') + outStr += printH(self.hy, outStr='\n hy:') + elif self.dim == 3: + outStr += '\n x0: {0:.2f}'.format(self.x0[0]) + outStr += '\n y0: {0:.2f}'.format(self.x0[1]) + outStr += '\n z0: {0:.2f}'.format(self.x0[2]) + outStr += '\n nCx: {0:d}'.format(self.nCx) + outStr += '\n nCy: {0:d}'.format(self.nCy) + outStr += '\n nCz: {0:d}'.format(self.nCz) + outStr += printH(self.hx, outStr='\n hx:') + outStr += printH(self.hy, outStr='\n hy:') + 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 diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 0205e972..02c23dee 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -100,11 +100,12 @@ except Exception, e: from InnerProducts import InnerProducts from TensorMesh import TensorMesh, BaseTensorMesh +from MeshIO import TreeMeshIO import time MAX_BITS = 20 -class TreeMesh(BaseTensorMesh, InnerProducts): +class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO): _meshType = 'TREE' diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index c43dfd13..eb5d13a1 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -102,408 +102,6 @@ def closestPoints(mesh, pts, gridLoc='CC'): return nodeInds -def readUBCTensorMesh(fileName): - """ - Read UBC GIF 3DTensor mesh and generate 3D Tensor mesh in simpegTD - - Input: - :param fileName, path to the UBC GIF mesh file - - Output: - :param SimPEG TensorMesh object - :return - """ - - # Interal function to read cell size lines for the UBC mesh files. - def readCellLine(line): - for seg in line.split(): - if '*' in seg: - st = seg - sp = seg.split('*') - re = np.array(sp[0],dtype=int)*(' ' + sp[1]) - line = line.replace(st,re.strip()) - return np.array(line.split(),dtype=float) - - # Read the file as line strings, remove lines with comment = ! - msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!') - - # Fist line is the size of the model - sizeM = np.array(msh[0].split(),dtype=float) - # Second line is the South-West-Top corner coordinates. - x0 = np.array(msh[1].split(),dtype=float) - # Read the cell sizes - h1 = readCellLine(msh[2]) - h2 = readCellLine(msh[3]) - h3temp = readCellLine(msh[4]) - h3 = h3temp[::-1] # Invert the indexing of the vector to start from the bottom. - # Adjust the reference point to the bottom south west corner - x0[2] = x0[2] - np.sum(h3) - # Make the mesh - from SimPEG import Mesh - tensMsh = Mesh.TensorMesh([h1,h2,h3],x0) - return tensMsh - -def readUBCTensorModel(fileName, mesh): - """ - Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg - - Input: - :param fileName, path to the UBC GIF mesh file to read - :param mesh, TensorMesh object, mesh that coresponds to the model - - Output: - :return numpy array, model with TensorMesh ordered - """ - f = open(fileName, 'r') - model = np.array(map(float, f.readlines())) - f.close() - model = np.reshape(model, (mesh.nCz, mesh.nCx, mesh.nCy), order = 'F') - model = model[::-1,:,:] - model = np.transpose(model, (1, 2, 0)) - model = mkvc(model) - - return model - -def writeUBCTensorMesh(fileName, mesh): - """ - Writes a SimPEG TensorMesh to a UBC-GIF format mesh file. - - :param str fileName: File to write to - :param simpeg.Mesh.TensorMesh mesh: The mesh - - """ - assert mesh.dim == 3 - s = '' - s += '%i %i %i\n' %tuple(mesh.vnC) - origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated. - origin.dtype = float - - 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[::-1]) - f = open(fileName, 'w') - f.write(s) - f.close() - -def writeUBCTensorModel(fileName, mesh, model): - """ - Writes a model associated with a SimPEG TensorMesh - to a UBC-GIF format model file. - - :param str fileName: File to write to - :param simpeg.Mesh.TensorMesh mesh: The mesh - :param numpy.ndarray model: The model - """ - - # Reshape model to a matrix - modelMat = mesh.r(model,'CC','CC','M') - # Transpose the axes - modelMatT = modelMat.transpose((2,0,1)) - # Flip z to positive down - modelMatTR = mkvc(modelMatT[::-1,:,:]) - - np.savetxt(fileName, modelMatTR.ravel()) - -def writeUBCocTreeFiles(fileName,mesh,modelDict=None): - ''' - Write UBC ocTree mesh and model files from a simpeg ocTree mesh and model. - - :param str fileName: File to write to - :param simpeg.Mesh.TreeMesh mesh: The mesh - :param dictionary modelDict: The models in a dictionary, where the keys is the name of the of the model file - - ''' - - # Calculate information to write in the file. - # Number of cells in the underlying mesh - nCunderMesh = np.array([h.size for h in mesh.h],dtype=np.int64) - # The top-south-west most corner of the mesh - tswCorn = mesh.x0 + np.array([0,0,np.sum(mesh.h[2])]) - # Smallest cell size - smallCell = np.array([h.min() for h in mesh.h]) - # Number of cells - nrCells = mesh.nC - - ## Extract iformation about the cells. - # cell pointers - cellPointers = np.array([c._pointer for c in mesh]) - # cell with - cellW = np.array([ mesh._levelWidth(i) for i in cellPointers[:,-1] ]) - # Need to shift the pointers to work with UBC indexing - # UBC Octree indexes always the top-left-close (top-south-west) corner first and orders the cells in z(top-down),x,y vs x,y,z(bottom-up). - # Shift index up by 1 - ubcCellPt = cellPointers[:,0:-1].copy() + np.array([1.,1.,1.]) - # Need reindex the z index to be from the top-left-close corner and to be from the global top. - ubcCellPt[:,2] = ( nCunderMesh[-1] + 2) - (ubcCellPt[:,2] + cellW) - - # Reorder the ubcCellPt - ubcReorder = np.argsort(ubcCellPt.view(','.join(3*['float'])),axis=0,order=['f2','f1','f0'])[:,0] - # Make a array with the pointers and the withs, that are order in the ubc ordering - indArr = np.concatenate((ubcCellPt[ubcReorder,:],cellW[ubcReorder].reshape((-1,1)) ),axis=1) - - ## Write the UBC octree mesh file - with open(fileName,'w') as mshOut: - mshOut.write('{:.0f} {:.0f} {:.0f}\n'.format(nCunderMesh[0],nCunderMesh[1],nCunderMesh[2])) - mshOut.write('{:.4f} {:.4f} {:.4f}\n'.format(tswCorn[0],tswCorn[1],tswCorn[2])) - mshOut.write('{:.3f} {:.3f} {:.3f}\n'.format(smallCell[0],smallCell[1],smallCell[2])) - mshOut.write('{:.0f} \n'.format(nrCells)) - np.savetxt(mshOut,indArr,fmt='%i') - - ## Print the models - # Assign the model('s) to the object - if modelDict is not None: - # indUBCvector = np.argsort(cX0[np.argsort(np.concatenate((cX0[:,0:2],cX0[:,2:3].max() - cX0[:,2:3]),axis=1).view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0]].view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0] - for item in modelDict.iteritems(): - # Save the data - np.savetxt(item[0],item[1][ubcReorder],fmt='%3.5e') - -def readUBCocTreeFiles(meshFile,modelFiles=None): - """ - Read UBC 3D OcTree mesh and/or modelFiles - - Input: - :param str meshFile: path to the UBC GIF OcTree mesh file to read - :param list of str modelFiles: list of paths modelFiles - - Output: - :return SimPEG.Mesh.TreeMesh mesh: The octree mesh - :return list of ndarray's: models as a list of numpy array's - """ - - ## Read the file lines - fileLines = np.genfromtxt(meshFile,dtype=str,delimiter='\n') - # Extract the data - nCunderMesh = np.array(fileLines[0].split(),dtype=float) - # I think this is the case? - if np.unique(nCunderMesh).size >1: - raise Exception('SimPEG TreeMeshes have the same number of cell in all directions') - tswCorn = np.array(fileLines[1].split(),dtype=float) - smallCell = np.array(fileLines[2].split(),dtype=float) - nrCells = np.array(fileLines[3].split(),dtype=float) - # Read the index array - indArr = np.genfromtxt(fileLines[4::],dtype=np.int) - - ## Calculate simpeg parameters - h1,h2,h3 = [np.ones(nr)*sz for nr,sz in zip(nCunderMesh,smallCell)] - x0 = tswCorn - np.array([0,0,np.sum(h3)]) - # Need to convert the index array to a points list that complies with SimPEG TreeMesh. - # Shift to start at 0 - simpegCellPt = indArr[:,0:-1].copy() - simpegCellPt[:,2] = ( nCunderMesh[-1] + 2) - (simpegCellPt[:,2] + indArr[:,3]) - # Need reindex the z index to be from the bottom-left-close corner and to be from the global bottom. - simpegCellPt = simpegCellPt - np.array([1.,1.,1.]) - - # Calculate the cell level - simpegLevel = np.log2(np.min(nCunderMesh)) - np.log2(indArr[:,3]) - # Make a pointer matrix - simpegPointers = np.concatenate((simpegCellPt,simpegLevel.reshape((-1,1))),axis=1) - - ## Make the tree mesh - from SimPEG.Mesh import TreeMesh - mesh = TreeMesh([h1,h2,h3],x0) - mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()]) - - # Figure out the reordering - simpegReorder = np.argsort(np.array([mesh._index(i) for i in simpegPointers.tolist()])) - # simpegReorder = np.argsort((np.array([[1,1,1,-1]])*simpegPointers).view(','.join(4*['float'])),axis=0,order=['f3','f2','f1','f0'])[:,0] - - if modelFiles is None: - return mesh - else: - modList = [] - for modFile in modelFiles: - modArr = np.loadtxt(modFile) - if len(modArr.shape) == 1: - modList.append(modArr[simpegReorder]) - else: - modList.append(modArr[simpegReorder,:]) - return mesh, modList - -def readVTRFile(fileName): - """ - Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model - - Input: - :param vtrFileName, path to the vtr model file to write to - - Output: - :return SimPEG TensorMesh object - :return SimPEG model dictionary - - """ - # Import - from vtk import vtkXMLRectilinearGridReader as vtrFileReader - from vtk.util.numpy_support import vtk_to_numpy - - # Read the file - vtrReader = vtrFileReader() - vtrReader.SetFileName(fileName) - vtrReader.Update() - vtrGrid = vtrReader.GetOutput() - # Sort information - hx = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetXCoordinates()))) - xR = vtk_to_numpy(vtrGrid.GetXCoordinates())[0] - hy = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetYCoordinates()))) - yR = vtk_to_numpy(vtrGrid.GetYCoordinates())[0] - zD = np.diff(vtk_to_numpy(vtrGrid.GetZCoordinates())) - # Check the direction of hz - if np.all(zD < 0): - hz = np.abs(zD[::-1]) - zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[-1] - else: - hz = np.abs(zD) - zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[0] - x0 = np.array([xR,yR,zR]) - - # Make the SimPEG object - from SimPEG import Mesh - tensMsh = Mesh.TensorMesh([hx,hy,hz],x0) - - # Grap the models - modelDict = {} - for i in np.arange(vtrGrid.GetCellData().GetNumberOfArrays()): - modelName = vtrGrid.GetCellData().GetArrayName(i) - if np.all(zD < 0): - modFlip = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i)) - tM = tensMsh.r(modFlip,'CC','CC','M') - modArr = tensMsh.r(tM[:,:,::-1],'CC','CC','V') - else: - modArr = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i)) - modelDict[modelName] = modArr - - # Return the data - return tensMsh, modelDict - -def writeVTRFile(fileName,mesh,model=None): - """ - Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model. - - Input: - :param str, path to the output vtk file - :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK - :param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells - - """ - # Import - from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter - from vtk.util.numpy_support import numpy_to_vtk - - # Deal with dimensionalities - if mesh.dim >= 1: - vX = mesh.vectorNx - xD = mesh.nNx - yD,zD = 1,1 - vY, vZ = np.array([0,0]) - if mesh.dim >= 2: - vY = mesh.vectorNy - yD = mesh.nNy - if mesh.dim == 3: - vZ = mesh.vectorNz - zD = mesh.nNz - # Use rectilinear VTK grid. - # Assign the spatial information. - vtkObj = rectGrid() - vtkObj.SetDimensions(xD,yD,zD) - vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1)) - vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1)) - vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1)) - - # Assign the model('s) to the object - if model is not None: - for item in model.iteritems(): - # Convert numpy array - vtkDoubleArr = numpy_to_vtk(item[1],deep=1) - vtkDoubleArr.SetName(item[0]) - vtkObj.GetCellData().AddArray(vtkDoubleArr) - # Set the active scalar - vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) - vtkObj.Update() - - - # Check the extension of the fileName - ext = os.path.splitext(fileName)[1] - if ext is '': - fileName = fileName + '.vtr' - elif ext not in '.vtr': - raise IOError('{:s} is an incorrect extension, has to be .vtr') - # Write the file. - vtrWriteFilter = rectWriter() - if float(VTK_VERSION.split('.')[0]) >=6: - vtrWriteFilter.SetInputData(vtkObj) - else: - vtuWriteFilter.SetInput(vtuObj) - vtrWriteFilter.SetFileName(fileName) - vtrWriteFilter.Update() - -def writeVTUFile(fileName,ocTreeMesh,modelDict=None): - ''' - Function to write a VTU file from a SimPEG TreeMesh and model. - ''' - from vtk import vtkXMLUnstructuredGridWriter as Writer, VTK_VERSION - from vtk.util.numpy_support import numpy_to_vtk - - # Make the object - vtuObj = simpegOcTree2vtuObj(ocTreeMesh,modelDict) - - # Make the writer - vtuWriteFilter = Writer() - if float(VTK_VERSION.split('.')[0]) >=6: - vtuWriteFilter.SetInputData(vtuObj) - else: - vtuWriteFilter.SetInput(vtuObj) - vtuWriteFilter.SetFileName(fileName) - # Write the file - vtuWriteFilter.Update() - -def simpegOcTree2vtuObj(simpegOcTreeMesh,modelDict=None): - ''' - Convert simpeg OcTree mesh and model to a VTK vtu object. - - ''' - import vtk - from vtk.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray - - if str(type(simpegOcTreeMesh)).split()[-1][1:-2] not in 'SimPEG.Mesh.TreeMesh.TreeMesh': - raise IOError('simpegOcTreeMesh is not a SimPEG TreeMesh.') - - # Make the data parts for the vtu object - # Points - try: - ptsMat = simpegOcTreeMesh._gridN + simpegOcTreeMesh.x0 - except: - simpegOcTreeMesh.number() - ptsMat = simpegOcTreeMesh._gridN + simpegOcTreeMesh.x0 - vtkPts = vtk.vtkPoints() - vtkPts.SetData(numpy_to_vtk(ptsMat,deep=True)) - # Cells - cellConn = np.array([c.nodes for c in simpegOcTreeMesh],dtype=np.int64) - - cellsMat = np.concatenate((np.ones((cellConn.shape[0],1),dtype=np.int64)*cellConn.shape[1],cellConn),axis=1).ravel() - cellsArr = vtk.vtkCellArray() - cellsArr.SetNumberOfCells(cellConn.shape[0]) - cellsArr.SetCells(cellConn.shape[0],numpy_to_vtkIdTypeArray(cellsMat,deep=True)) - - # Make the object - vtuObj = vtk.vtkUnstructuredGrid() - vtuObj.SetPoints(vtkPts) - vtuObj.SetCells(vtk.VTK_VOXEL,cellsArr) - # Add the level of refinement as a cell array - cellSides = np.array([np.array(vtuObj.GetCell(i).GetBounds()).reshape((3,2)).dot(np.array([-1, 1])) for i in np.arange(vtuObj.GetNumberOfCells())]) - uniqueLevel, indLevel = np.unique(np.prod(cellSides,axis=1),return_inverse=True) - refineLevelArr = numpy_to_vtk(indLevel.max() - indLevel,deep=1) - refineLevelArr.SetName('octreeLevel') - vtuObj.GetCellData().AddArray(refineLevelArr) - # Assign the model('s) to the object - if modelDict is not None: - for item in modelDict.iteritems(): - # Convert numpy array - vtkDoubleArr = numpy_to_vtk(item[1],deep=1) - vtkDoubleArr.SetName(item[0]) - vtuObj.GetCellData().AddArray(vtkDoubleArr) - - return vtuObj - def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'): """ Extracts Core Mesh from Global mesh diff --git a/tests/mesh/test_MeshIO.py b/tests/mesh/test_MeshIO.py index 52e5740d..e8dc0748 100644 --- a/tests/mesh/test_MeshIO.py +++ b/tests/mesh/test_MeshIO.py @@ -4,11 +4,65 @@ import SimPEG as simpeg from SimPEG.Mesh import TensorMesh, TreeMesh -class TestOcTreeIO(unittest.TestCase): +class TestTensorMeshIO(unittest.TestCase): def setUp(self): h = np.ones(16) - mesh = simpeg.Mesh.TreeMesh([h,2*h,3*h]) + 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]) @@ -19,9 +73,10 @@ class TestOcTreeIO(unittest.TestCase): mesh = self.mesh # Make a vector vec = np.arange(mesh.nC) - # Write aand read - simpeg.Utils.meshutils.writeUBCocTreeFiles('temp.msh',mesh,{'arange.txt':vec}) - meshUBC, vecUBC = simpeg.Utils.meshutils.readUBCocTreeFiles('temp.msh',['arange.txt']) + # 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__() @@ -35,7 +90,7 @@ class TestOcTreeIO(unittest.TestCase): def test_VTUfiles(self): mesh = self.mesh vec = np.arange(mesh.nC) - simpeg.Utils.meshutils.writeVTUFile('temp.vtu',mesh,{'arange':vec}) + mesh.writeVTK('temp.vtu',{'arange':vec}) print 'Writing of VTU files is working' os.remove('temp.vtu') From b0135082707961b27f5a45cf3ca872b2ad0ed019 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 28 Jan 2016 22:33:04 -0800 Subject: [PATCH 30/30] Fixes #215 --- SimPEG/EM/TDEM/BaseTDEM.py | 10 +++++++++- SimPEG/Problem.py | 12 ++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/SimPEG/EM/TDEM/BaseTDEM.py b/SimPEG/EM/TDEM/BaseTDEM.py index 4d4b39a7..2efb10ec 100644 --- a/SimPEG/EM/TDEM/BaseTDEM.py +++ b/SimPEG/EM/TDEM/BaseTDEM.py @@ -37,13 +37,21 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): _FieldsForward_pair = FieldsTDEM #: used for the forward calculation only + waveformType = "STEPOFF" + current = None + + def currentwaveform(self, wave): + self._timeSteps = np.diff(wave[:,0]) + self.current = wave[:,1] + self.waveformType = "GENERAL" + def fields(self, m): if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50) self.curModel = m # Create a fields storage object F = self._FieldsForward_pair(self.mesh, self.survey) for src in self.survey.srcList: - # Set the initial conditions + # Set the initial conditions F[src,:,0] = src.getInitialFields(self.mesh) F = self.forward(m, self.getRHS, F=F) if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index be29ff01..f24b57de 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -32,8 +32,8 @@ class BaseProblem(object): val._assertMatchesPair(self.mapPair) self._mapping = val else: - self._mapping = self.PropMap(val) - + self._mapping = self.PropMap(val) + def __init__(self, mesh, mapping=None, **kwargs): Utils.setKwargs(self, **kwargs) assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object." @@ -158,9 +158,6 @@ class BaseProblem(object): class BaseTimeProblem(BaseProblem): """Sets up that basic needs of a time domain problem.""" - - waveformType = "STEPOFF" - current = None @property def timeSteps(self): @@ -187,11 +184,6 @@ class BaseTimeProblem(BaseProblem): self._timeSteps = Utils.meshTensor(value) del self.timeMesh - def currentwaveform(self, wave): - self._timeSteps = np.diff(wave[:,0]) - self.current = wave[:,1] - self.waveformType = "GENERAL" - @property def nT(self): "Number of time steps."