Fixed bugs in meshutils, writing VTR files

Allow fileName to be None, which will output the VTKobject without save a file.

SimPEG/DataMisfit.py

@@ -14,6 +14,34 @@ class BaseDataMisfit(object):
     debug   = False  #: Print debugging information
     counter = None   #: Set this to a SimPEG.Utils.Counter() if you want to count things
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     def __init__(self, survey, **kwargs):
         assert survey.ispaired, 'The survey must be paired to a problem.'
         if isinstance(survey, Survey.BaseSurvey):

SimPEG/Directives.py

@@ -8,6 +8,34 @@ class InversionDirective(object):
     def __init__(self, **kwargs):
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def inversion(self):
         """This is the inversion of the InversionDirective instance."""
@@ -207,6 +235,76 @@ class SaveOutputEveryIteration(_SaveEveryIteration):
         f.close()
 
 
+class SaveOutputDictEveryIteration(_SaveEveryIteration):
+    """SaveOutputDictEveryIteration"""
+
+    def initialize(self):
+        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.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) )
+        phi_mx = 0.5 * mx.dot(mx)
+        if self.prob.mesh.dim==2:
+            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 - 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 SaveOutputDictEveryIteration(_SaveEveryIteration):
+    """SaveOutputDictEveryIteration
+
+    A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
+
+
+
+    """
+
+    def initialize(self):
+        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.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) )
+        phi_mx = 0.5 * mx.dot(mx)
+        if self.prob.mesh.dim==2:
+            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 - 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,dpred=self.invProb.dpred)
+
+
 
 
 # class UpdateReferenceModel(Parameter):

SimPEG/Fields.py

@@ -12,6 +12,34 @@ class Fields(object):
     aliasFields = None  #: Aliased fields, a dict with [alias, location, function], e.g. {"b":["e","F",lambda(F,e,ind)]}
     dtype = float       #: dtype is the type of the storage matrix. This can be a dictionary.
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     def __init__(self, mesh, survey, **kwargs):
         self.survey = survey
         self.mesh = mesh

SimPEG/Maps.py

@@ -17,6 +17,30 @@ class IdentityMap(object):
         Utils.setKwargs(self, **kwargs)
         self.mesh = mesh
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+
     @property
     def nP(self):
         """
@@ -289,7 +313,7 @@ class FullMap(IdentityMap):
     """
     FullMap
 
-    Given a scalar, the FullMap maps the value to the 
+    Given a scalar, the FullMap maps the value to the
     full model space.
     """
 
@@ -314,8 +338,8 @@ class FullMap(IdentityMap):
             :rtype: numpy.array
             :return: derivative of transformed model
         """
-        return np.ones([self.mesh.nC,1])     
-      
+        return np.ones([self.mesh.nC,1])
+
 
 class Vertical1DMap(IdentityMap):
     """Vertical1DMap

SimPEG/Mesh/InnerProducts.py

@@ -33,7 +33,7 @@ class InnerProducts(object):
         return self._getInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat, doFast=doFast)
 
     def _getInnerProduct(self, projType, prop=None, invProp=False, invMat=False, doFast=True):
-        """
+        """r
             :param str projType: 'F' for faces 'E' for edges
             :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param bool invProp: inverts the material property

SimPEG/Optimization.py

@@ -115,6 +115,34 @@ class Minimize(object):
 
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def callback(self):
         return getattr(self, '_callback', None)

SimPEG/Problem.py

@@ -22,6 +22,34 @@ class BaseProblem(object):
 
     PropMap = None    #: A SimPEG PropertyMap class.
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def mapping(self):
         "A SimPEG.Map instance or a property map is PropMap is not None"
@@ -32,8 +60,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."

SimPEG/PropMaps.py

@@ -13,6 +13,34 @@ class Property(object):
         self.doc = doc
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def propertyLink(self):
         "Can be something like: ('sigma', Maps.ReciprocalMap)"
@@ -118,6 +146,36 @@ class PropModel(object):
         self.vector  = vector
         assert len(self.vector) == self.nP
 
+    # Pickleing support methods
+    # def __reduce__(self):
+    #     return (dict,{self.propMap,self.vector})
+
+    # def __getstate__(self):
+    #     '''
+    #     Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+    #     Used when doing:
+    #         pickle.dump(pickleFile,object)
+    #     '''
+    #     self.__class__ = ProbModel
+    #     odict = {}
+    #     odict['vec'] = self.__dict__['vector']
+    #     odict['pMap'] = self.__dict__['propMap']
+    #     # Return the dict
+    #     return odict
+
+    # def __setstate__(self,odict):
+    #     '''
+    #     Function that sets a pickle dictionary in to an object.
+
+    #     Used when doing:
+    #         object = pickle.load(pickleFile)
+    #     '''
+    #     # Update the dict
+    #     # Re-hook the methods to the object
+    #     self.propMap = odict['prMap']
+    #     self.vector = odict['vec']
+
     @property
     def nP(self):
         inds = []
@@ -207,6 +265,24 @@ class PropMap(object):
         else:
             raise Exception('mappings must be a dict, a mapping, or a list of tuples.')
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        pass
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        pass
 
     def setup(self, maps, slices=None):
         """
@@ -239,7 +315,7 @@ class PropMap(object):
             setattr(self, '%sMap'%name, mapping)
             setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP)))
             nP += mapping.nP
-        self.nP = nP 
+        self.nP = nP
 
     @property
     def defaultInvProp(self):

SimPEG/Regularization.py

@@ -27,6 +27,34 @@ class BaseRegularization(object):
         self.mapping = mapping or Maps.IdentityMap(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def parent(self):
         """This is the parent of the regularization."""
@@ -311,7 +339,7 @@ class Tikhonov(BaseRegularization):
         if self.smoothModel == True:
             mD1 = self.mapping.deriv(m)
             mD2 = self.mapping.deriv(m - self.mref)
-            r1 = self.Wsmooth * ( self.mapping * (m)) 
+            r1 = self.Wsmooth * ( self.mapping * (m))
             r2 = self.Ws * ( self.mapping * (m - self.mref) )
             out1 = mD1.T * ( self.Wsmooth.T * r1 )
             out2 = mD2.T * ( self.Ws.T * r2 )

SimPEG/Survey.py

@@ -19,6 +19,35 @@ class BaseRx(object):
         self._Ps = {}
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
+
     @property
     def rxType(self):
         """Receiver Type"""
@@ -129,6 +158,33 @@ class BaseSrc(object):
         self.rxList = rxList
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
 
     @property
     def nD(self):
@@ -153,6 +209,25 @@ class Data(object):
         if v is not None:
             self.fromvec(v)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        pass
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        pass
+
     def _ensureCorrectKey(self, key):
         if type(key) is tuple:
             if len(key) is not 2:
@@ -210,11 +285,39 @@ class BaseSurvey(object):
     mtrue = None     #: True model, if data is synthetic
 
     counter = None   #: A SimPEG.Utils.Counter object
+    srcPair = BaseSrc  #: Source Pair
 
     def __init__(self, **kwargs):
         Utils.setKwargs(self, **kwargs)
 
-    srcPair = BaseSrc  #: Source Pair
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
 
     @property
     def srcList(self):

SimPEG/Utils/ModelBuilder.py

@@ -170,16 +170,58 @@ def scalarConductivity(ccMesh,pFunction):
 
     return mkvc(sigma)
 
+def layeredModel(ccMesh, layerTops, layerValues):
+    """
+        Define a layered model from layerTops (z-positive up)
+
+        :param numpy.array ccMesh: cell-centered mesh
+        :param numpy.array layerTops: z-locations of the tops of each layer
+        :param numpy.array layerValue: values of the property to assign for each layer (starting at the top) 
+        :rtype: numpy.array
+        :return: M, layered model on the mesh 
+    """
+
+    descending = np.linalg.norm(sorted(layerTops, reverse=True) - layerTops) < 1e-20
+
+    # TODO: put an error check to make sure that there is an ordering... needs to work with inf elts 
+    # assert ascending or descending, "Layers must be listed in either ascending or descending order"
+
+    # start from bottom up
+    if not descending:
+        zprop = np.hstack([mkvc(layerTops,2),mkvc(layerValues,2)])
+        zprop.sort(axis=0)
+        layerTops, layerValues  = zprop[::-1,0], zprop[::-1,1]
+
+    # put in vector form
+    layerTops, layerValues = mkvc(layerTops), mkvc(layerValues)
+
+    # initialize with bottom layer
+    dim = ccMesh.shape[1]
+    if dim == 3:
+        z = ccMesh[:,2]
+    elif dim == 2:
+        z = ccMesh[:,1]
+    elif dim == 1:
+        z = ccMesh[:,0]
+
+    model = np.zeros(ccMesh.shape[0])
+
+    for i, top in enumerate(layerTops):
+        zind = z <= top 
+        model[zind] = layerValues[i]
+
+    return model 
+
 
 
 def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
     """
-        Create a random model by convolving a kernal with a
+        Create a random model by convolving a kernel with a
         uniformly distributed model.
 
         :param int,tuple shape: shape of the model.
         :param int seed: pick which model to produce, prints the seed if you don't choose.
-        :param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernal that is used.
+        :param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernel that is used.
         :param int its: number of smoothing iterations
         :param list bounds: bounds on the model, len(list) == 2
         :rtype: numpy.ndarray

SimPEG/Utils/meshutils.py

@@ -149,7 +149,7 @@ def readUBCTensorModel(fileName, mesh):
 
         Input:
         :param fileName, path to the UBC GIF mesh file to read
-        :param mesh, TensorMesh object, mesh that coresponds to the model 
+        :param mesh, TensorMesh object, mesh that coresponds to the model
 
         Output:
         :return numpy array, model with TensorMesh ordered
@@ -170,7 +170,7 @@ def writeUBCTensorMesh(fileName, mesh):
 
         :param str fileName: File to write to
         :param simpeg.Mesh.TensorMesh mesh: The mesh
-        
+
     """
     assert mesh.dim == 3
     s = ''
@@ -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
@@ -216,7 +215,7 @@ def readVTRFile(fileName):
         Output:
         :return SimPEG TensorMesh object
         :return SimPEG model dictionary
-        
+
     """
     # Import
     from vtk import vtkXMLRectilinearGridReader as vtrFileReader
@@ -296,84 +295,87 @@ 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])
-    vtkObj.Update()
+    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])
 
 
     # 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()
-    vtrWriteFilter.SetInput(vtkObj)
-    vtrWriteFilter.SetFileName(fileName)
-    vtrWriteFilter.Update()
+    if fileName is not None:
+        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()
+        vtrWriteFilter.SetInputData(vtkObj)
+        vtrWriteFilter.SetFileName(fileName)
+        vtrWriteFilter.Update()
+    else:
+        return vtkObj
 
 def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
     """
         Extracts Core Mesh from Global mesh
         xyzlim: 2D array [ndim x 2]
         mesh: SimPEG mesh
-        This function ouputs: 
+        This function ouputs:
             - actind: corresponding boolean index from global to core
-            - meshcore: core SimPEG mesh  
+            - meshcore: core SimPEG mesh
         Warning: 1D and 2D has not been tested
     """
     from SimPEG import Mesh
     if mesh.dim ==1:
         xyzlim = xyzlim.flatten()
         xmin, xmax = xyzlim[0], xyzlim[1]
-        
-        xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)        
-        
+
+        xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)
+
         xc = mesh.vectorCCx[xind]
 
         hx = mesh.hx[xind]
-        
+
         x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
-        
+
         meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
-        
+
         actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax)
-        
+
     elif mesh.dim ==2:
         xmin, xmax = xyzlim[0,0], xyzlim[0,1]
         ymin, ymax = xyzlim[1,0], xyzlim[1,1]
 
         yind = np.logical_and(mesh.vectorCCy>ymin, mesh.vectorCCy<ymax)
-        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)        
+        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)
 
         xc = mesh.vectorCCx[xind]
         yc = mesh.vectorCCy[yind]
 
         hx = mesh.hx[xind]
         hy = mesh.hy[yind]
-        
+
         x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
-        
+
         meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
-        
+
         actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
                & (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
-        
+
     elif mesh.dim==3:
         xmin, xmax = xyzlim[0,0], xyzlim[0,1]
         ymin, ymax = xyzlim[1,0], xyzlim[1,1]
         zmin, zmax = xyzlim[2,0], xyzlim[2,1]
-        
+
         xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)
         yind = np.logical_and(mesh.vectorCCy>ymin, mesh.vectorCCy<ymax)
-        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)        
+        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)
 
         xc = mesh.vectorCCx[xind]
         yc = mesh.vectorCCy[yind]
@@ -382,19 +384,19 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
         hx = mesh.hx[xind]
         hy = mesh.hy[yind]
         hz = mesh.hz[zind]
-        
+
         x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5, zc[0]-hz[0]*0.5]
-        
+
         meshCore = Mesh.TensorMesh([hx, hy, hz] ,x0=x0)
-        
+
         actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
                & (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
                & (mesh.gridCC[:,2]>zmin) & (mesh.gridCC[:,2]<zmax)
-                
+
     else:
         raise(Exception("Not implemented!"))
-    
-    
+
+
     return actind, meshCore