diff --git a/.bumpversion.cfg b/.bumpversion.cfg
index 189947df..f739d080 100644
--- a/.bumpversion.cfg
+++ b/.bumpversion.cfg
@@ -1,4 +1,4 @@
 [bumpversion]
-current_version = 0.1.1
+current_version = 0.1.3
 files = setup.py SimPEG/__init__.py docs/conf.py
 
diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py
index e84a7d39..3ec9e96e 100644
--- a/SimPEG/Directives.py
+++ b/SimPEG/Directives.py
@@ -144,23 +144,69 @@ class BetaSchedule(InversionDirective):
             if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
             self.invProb.beta /= self.coolingFactor
 
-
-
-class SaveModelEveryIteration(InversionDirective):
-    """SaveModelEveryIteration"""
+class TargetMisfit(InversionDirective):
 
     @property
-    def modelName(self):
-        if getattr(self, '_modelName', None) is None:
-            from datetime import datetime
-            self._modelName = 'inversionModel-%s'%datetime.now().strftime('%Y-%m-%d')
-        return self._modelName
-    @modelName.setter
-    def modelName(self, value):
-        self._modelName = value
+    def target(self):
+        if getattr(self, '_target', None) is None:
+            self._target = self.survey.nD
+        return self._target
+    @target.setter
+    def target(self, val):
+        self._target = val
 
     def endIter(self):
-        np.save('%03d-%s' % (self.opt.iter, self.modelName), self.opt.xc)
+        if self.invProb.phi_d < self.target:
+            self.opt.stopNextIteration = True
+
+
+
+class _SaveEveryIteration(InversionDirective):
+    @property
+    def name(self):
+        if getattr(self, '_name', None) is None:
+            self._name = 'InversionModel'
+        return self._name
+    @name.setter
+    def name(self, value):
+        self._name = value
+
+    @property
+    def fileName(self):
+        if getattr(self, '_fileName', None) is None:
+            from datetime import datetime
+            self._fileName = '%s-%s'%(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M'))
+        return self._fileName
+    @fileName.setter
+    def fileName(self, value):
+        self._fileName = value
+
+
+class SaveModelEveryIteration(_SaveEveryIteration):
+    """SaveModelEveryIteration"""
+
+    def initialize(self):
+        print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName
+
+    def endIter(self):
+        np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
+
+
+class SaveOutputEveryIteration(_SaveEveryIteration):
+    """SaveModelEveryIteration"""
+
+    def initialize(self):
+        print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName
+        f = open(self.fileName+'.txt', 'w')
+        f.write("  #     beta     phi_d     phi_m       f\n")
+        f.close()
+
+    def endIter(self):
+        f = open(self.fileName+'.txt', 'a')
+        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 UpdateReferenceModel(Parameter):
diff --git a/SimPEG/Examples/DCfwd.py b/SimPEG/Examples/DCfwd.py
new file mode 100644
index 00000000..385babf6
--- /dev/null
+++ b/SimPEG/Examples/DCfwd.py
@@ -0,0 +1,71 @@
+from SimPEG import Mesh, Utils, np, SolverLU
+import matplotlib.pyplot as plt
+import matplotlib
+from matplotlib.mlab import griddata
+
+## 2D DC forward modeling example with Tensor and Curvilinear Meshes
+
+# Step1: Generate Tensor and Curvilinear Mesh
+sz = [40,40]
+# Tensor Mesh
+tM = Mesh.TensorMesh(sz)
+# Curvilinear Mesh
+rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
+
+# Step2: Direct Current (DC) operator
+def DCfun(mesh, pts):
+    D = mesh.faceDiv
+    G = D.T
+    sigma = 1e-2*np.ones(mesh.nC)
+    Msigi = mesh.getFaceInnerProduct(1./sigma)
+    MsigI = Utils.sdInv(Msigi)
+    A = D*MsigI*G
+    A[-1,-1] /= mesh.vol[-1] # Remove null space
+    rhs = np.zeros(mesh.nC)
+    txind = Utils.meshutils.closestPoints(mesh, pts)
+    rhs[txind] = np.r_[1,-1]
+    return A, rhs
+
+pts = np.vstack((np.r_[0.25, 0.5], np.r_[0.75, 0.5]))
+
+#Step3: Solve DC problem (LU solver)
+AtM, rhstM = DCfun(tM, pts)
+AinvtM = SolverLU(AtM)
+phitM = AinvtM*rhstM
+
+ArM, rhsrM = DCfun(rM, pts)
+AinvrM = SolverLU(ArM)
+phirM = AinvrM*rhsrM
+
+#Step4: Making Figure
+fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
+label = ["(a)", "(b)"]
+opts = {}
+vmin, vmax = phitM.min(), phitM.max()
+dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True)
+
+#TODO: At the moment Curvilinear Mesh do not have plotimage 
+
+Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F')
+Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F')
+PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear')
+axes[1].contourf(Xi, Yi, PHIrM, 100, vmin=vmin, vmax=vmax)
+
+cb = plt.colorbar(dat[0], ax=axes[0]); cb.set_label("Voltage (V)")
+cb = plt.colorbar(dat[0], ax=axes[1]); cb.set_label("Voltage (V)")
+
+tM.plotGrid(ax=axes[0], **opts)
+axes[0].set_title('TensorMesh')
+rM.plotGrid(ax=axes[1], **opts)
+axes[1].set_title('CurvilinearMesh')
+for i in range(2):
+    axes[i].set_xlim(0.025, 0.975)
+    axes[i].set_ylim(0.025, 0.975)
+    axes[i].text(0., 1.0, label[i], fontsize=20)
+    if i==0:        
+        axes[i].set_ylabel("y")
+    else:
+        axes[i].set_ylabel(" ")
+    axes[i].set_xlabel("x")
+
+plt.show()
diff --git a/SimPEG/Fields.py b/SimPEG/Fields.py
new file mode 100644
index 00000000..9baa32f8
--- /dev/null
+++ b/SimPEG/Fields.py
@@ -0,0 +1,261 @@
+import Utils, numpy as np, scipy.sparse as sp
+
+class Fields(object):
+    """Fancy Field Storage
+
+        u[:,'phi'] = phi
+        print u[src0,'phi']
+
+    """
+
+    knownFields = None  #: Known fields,   a dict with locations,         e.g. {"e": "E", "phi": "CC"}
+    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.
+
+    def __init__(self, mesh, survey, **kwargs):
+        self.survey = survey
+        self.mesh = mesh
+        Utils.setKwargs(self, **kwargs)
+        self._fields = {}
+
+        if self.knownFields is None:
+            raise Exception('knownFields cannot be set to None')
+        if self.aliasFields is None:
+            self.aliasFields = {}
+
+        allFields = [k for k in self.knownFields] + [a for a in self.aliasFields]
+        assert len(allFields) == len(set(allFields)), 'Aliased fields and Known Fields have overlapping definitions.'
+        self.startup()
+
+    def startup(self):
+        pass
+
+    @property
+    def approxSize(self):
+        """The approximate cost to storing all of the known fields."""
+        sz = 0.0
+        for f in self.knownFields:
+            loc =self.knownFields[f]
+            sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2)
+        return "%e MB"%sz
+
+    def _storageShape(self, loc):
+        nSrc = self.survey.nSrc
+
+        nP = {'CC': self.mesh.nC,
+              'N':  self.mesh.nN,
+              'F':  self.mesh.nF,
+              'E':  self.mesh.nE}[loc]
+
+        return (nP, nSrc)
+
+    def _initStore(self, name):
+        if name in self._fields:
+            return self._fields[name]
+
+        assert name in self.knownFields, 'field name is not known.'
+
+        loc = self.knownFields[name]
+
+        if type(self.dtype) is dict:
+            dtype = self.dtype[name]
+        else:
+            dtype = self.dtype
+        field = np.zeros(self._storageShape(loc), dtype=dtype)
+
+        self._fields[name] = field
+
+        return field
+
+    def _srcIndex(self, srcTestList):
+        if type(srcTestList) is slice:
+            ind = srcTestList
+        else:
+            ind = self.survey.getSourceIndex(srcTestList)
+        return ind
+
+    def _nameIndex(self, name, accessType):
+
+        if type(name) is slice:
+            assert name == slice(None,None,None), 'Fancy field name slicing is not supported... yet.'
+            name = None
+
+        if name is None:
+            return
+        if accessType=='set' and name not in self.knownFields:
+            if name in self.aliasFields:
+                raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name)
+            else:
+                raise KeyError('Invalid field name (%s) for setter'%name)
+
+        elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields):
+            raise KeyError('Invalid field name (%s) for getter'%name)
+        return name
+
+    def _indexAndNameFromKey(self, key, accessType):
+        if type(key) is not tuple:
+            key = (key,)
+        if len(key) == 1:
+            key += (None,)
+
+        assert len(key) == 2, 'must be [Src, fieldName]'
+
+        srcTestList, name = key
+        name = self._nameIndex(name, accessType)
+        ind = self._srcIndex(srcTestList)
+        return ind, name
+
+    def __setitem__(self, key, value):
+        ind, name = self._indexAndNameFromKey(key, 'set')
+        if name is None:
+            freq = key
+            assert type(value) is dict, 'New fields must be a dictionary, if field is not specified.'
+            newFields = value
+        elif name in self.knownFields:
+            newFields = {name: value}
+        else:
+            raise Exception('Unknown setter')
+
+        for name in newFields:
+            field = self._initStore(name)
+            self._setField(field, newFields[name], name, ind)
+
+    def __getitem__(self, key):
+        ind, name = self._indexAndNameFromKey(key, 'get')
+        if name is None:
+            out = {}
+            for name in self._fields:
+                out[name] = self._getField(name, ind)
+            return out
+        return self._getField(name, ind)
+
+    def _setField(self, field, val, name, ind):
+        if isinstance(val, np.ndarray) and (field.shape[0] == field.size or val.ndim == 1):
+            val = Utils.mkvc(val,2)
+        field[:,ind] = val
+
+    def _getField(self, name, ind):
+        if name in self._fields:
+            out = self._fields[name][:,ind]
+        else:
+            # Aliased fields
+            alias, loc, func = self.aliasFields[name]
+
+            srcII = np.array(self.survey.srcList)[ind]
+            srcII = srcII.tolist()
+
+            if type(func) is str:
+                assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
+                func = getattr(self, func)
+            out = func(self._fields[alias][:,ind], srcII)
+        if out.shape[0] == out.size or out.ndim == 1:
+            out = Utils.mkvc(out,2)
+        return out
+
+    def __contains__(self, other):
+        if other in self.aliasFields:
+            other = self.aliasFields[other][0]
+        return self._fields.__contains__(other)
+
+
+class TimeFields(Fields):
+    """Fancy Field Storage for time domain problems
+
+        u[:,'phi', timeInd] = phi
+        print u[src0,'phi']
+
+    """
+
+    def _storageShape(self, loc):
+        nP = {'CC': self.mesh.nC,
+              'N':  self.mesh.nN,
+              'F':  self.mesh.nF,
+              'E':  self.mesh.nE}[loc]
+        nSrc = self.survey.nSrc
+        nT = self.survey.prob.nT + 1
+        return (nP, nSrc, nT)
+
+    def _indexAndNameFromKey(self, key, accessType):
+        if type(key) is not tuple:
+            key = (key,)
+        if len(key) == 1:
+            key += (None,)
+        if len(key) == 2:
+            key += (slice(None,None,None),)
+
+        assert len(key) == 3, 'must be [Src, fieldName, times]'
+
+        srcTestList, name, timeInd = key
+
+        name = self._nameIndex(name, accessType)
+        srcInd = self._srcIndex(srcTestList)
+
+        return (srcInd, timeInd), name
+
+    def _correctShape(self, name, ind, deflate=False):
+        srcInd, timeInd = ind
+        if name in self.knownFields:
+            loc = self.knownFields[name]
+        else:
+            loc = self.aliasFields[name][1]
+        nP, total_nSrc, total_nT = self._storageShape(loc)
+        nSrc = np.ones(total_nSrc, dtype=bool)[srcInd].sum()
+        nT  = np.ones(total_nT, dtype=bool)[timeInd].sum()
+        shape = nP, nSrc, nT
+        if deflate:
+             shape = tuple([s for s in shape if s > 1])
+        if len(shape) == 1:
+            shape = shape + (1,)
+        return shape
+
+    def _setField(self, field, val, name, ind):
+        srcInd, timeInd = ind
+        shape = self._correctShape(name, ind)
+        if Utils.isScalar(val):
+            field[:,srcInd,timeInd] = val
+            return
+        if val.size != np.array(shape).prod():
+            raise ValueError('Incorrect size for data.')
+        correctShape = field[:,srcInd,timeInd].shape
+        field[:,srcInd,timeInd] = val.reshape(correctShape, order='F')
+
+    def _getField(self, name, ind):
+        srcInd, timeInd = ind
+
+        if name in self._fields:
+            out = self._fields[name][:,srcInd,timeInd]
+        else:
+            # Aliased fields
+            alias, loc, func = self.aliasFields[name]
+            if type(func) is str:
+                assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
+                func = getattr(self, func)
+            pointerFields = self._fields[alias][:,srcInd,timeInd]
+            pointerShape = self._correctShape(alias, ind)
+            pointerFields = pointerFields.reshape(pointerShape, order='F')
+
+            timeII = np.arange(self.survey.prob.nT + 1)[timeInd]
+            srcII  = np.array(self.survey.srcList)[srcInd]
+            srcII  = srcII.tolist()
+
+            if timeII.size == 1:
+                pointerShapeDeflated = self._correctShape(alias, ind, deflate=True)
+                pointerFields = pointerFields.reshape(pointerShapeDeflated, order='F')
+                out = func(pointerFields, srcII, timeII)
+            else: #loop over the time steps
+                nT = pointerShape[2]
+                out = range(nT)
+                for i, TIND_i in enumerate(timeII):
+                    fieldI = pointerFields[:,:,i]
+                    if fieldI.shape[0] == fieldI.size:
+                        fieldI = Utils.mkvc(fieldI, 2)
+                    out[i] = func(fieldI, srcII, TIND_i)
+                    if out[i].ndim == 1:
+                        out[i] = out[i][:,np.newaxis,np.newaxis]
+                    elif out[i].ndim == 2:
+                        out[i] = out[i][:,:,np.newaxis]
+                out = np.concatenate(out, axis=2)
+
+        shape = self._correctShape(name, ind, deflate=True)
+        return out.reshape(shape, order='F')
+
diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py
index 0b47de66..74a92efd 100644
--- a/SimPEG/Maps.py
+++ b/SimPEG/Maps.py
@@ -1,5 +1,6 @@
 import Utils, numpy as np, scipy.sparse as sp
 from Tests import checkDerivative
+from PropMaps import PropMap, Property
 
 
 class IdentityMap(object):
@@ -22,6 +23,8 @@ class IdentityMap(object):
             :rtype: int
             :return: number of parameters in the model
         """
+        if self.mesh is None:
+            return '*'
         return self.mesh.nC
 
     @property
@@ -32,8 +35,11 @@ class IdentityMap(object):
             :rtype: (int,int)
             :return: shape of the operator as a tuple
         """
+        if self.mesh is None:
+            return ('*', self.nP)
         return (self.mesh.nC, self.nP)
 
+
     def _transform(self, m):
         """
             Changes the model into the physical property.
@@ -98,17 +104,17 @@ class IdentityMap(object):
 
     def __mul__(self, val):
         if isinstance(val, IdentityMap):
-            if not self.shape[1] == val.shape[0]:
+            if not (self.shape[1] == '*' or val.shape[0] == '*') and not self.shape[1] == val.shape[0]:
                 raise ValueError('Dimension mismatch in %s and %s.' % (str(self), str(val)))
             return ComboMap([self, val])
         elif isinstance(val, np.ndarray):
-            if not self.shape[1] == val.shape[0]:
+            if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]:
                 raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape)))
             return self._transform(val)
         raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!')
 
     def __str__(self):
-        return "%s(%d,%d)" % (self.__class__.__name__, self.shape[0], self.shape[1])
+        return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
 
 class ComboMap(IdentityMap):
     """Combination of various maps."""
@@ -119,10 +125,10 @@ class ComboMap(IdentityMap):
         self.maps = []
         for ii, m in enumerate(maps):
             assert isinstance(m, IdentityMap), 'Unrecognized data type, inherit from an IdentityMap or ComboMap!'
-            if ii > 0 and not self.shape[1] == m.shape[0]:
+            if ii > 0 and not (self.shape[1] == '*' or m.shape[0] == '*') and not self.shape[1] == m.shape[0]:
                 prev = self.maps[-1]
-                errArgs = (prev.__name__, prev.shape[0], prev.shape[1], m.__name__, m.shape[0], m.shape[1])
-                raise ValueError('Dimension mismatch in map[%s] (%i, %i) and map[%s] (%i, %i).' % errArgs)
+                errArgs = (prev.__class__.__name__, prev.shape[0], prev.shape[1], m.__class__.__name__, m.shape[0], m.shape[1])
+                raise ValueError('Dimension mismatch in map[%s] (%s, %s) and map[%s] (%s, %s).' % errArgs)
 
             if isinstance(m, ComboMap):
                 self.maps += m.maps
@@ -155,7 +161,7 @@ class ComboMap(IdentityMap):
         return deriv
 
     def __str__(self):
-        return 'ComboMap[%s]%s' % (' * '.join([m.__str__() for m in self.maps]), str(self.shape))
+        return 'ComboMap[%s](%s,%s)' % (' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1])
 
 
 class ExpMap(IdentityMap):
@@ -220,6 +226,26 @@ class ExpMap(IdentityMap):
         """
         return Utils.sdiag(np.exp(Utils.mkvc(m)))
 
+class ReciprocalMap(IdentityMap):
+    """
+        Reciprocal mapping. For example, electrical resistivity and conductivity.
+
+        .. math::
+
+            \\rho = \\frac{1}{\sigma}
+
+    """
+    def _transform(self, m):
+        return 1.0 / Utils.mkvc(m)
+
+    def inverse(self, D):
+        return 1.0 / Utils.mkvc(m)
+
+    def deriv(self, m):
+        # TODO: if this is a tensor, you might have a problem.
+        return Utils.sdiag( - Utils.mkvc(m)**(-2) )
+
+
 
 class LogMap(IdentityMap):
     """
diff --git a/SimPEG/Mesh/LogicallyRectMesh.py b/SimPEG/Mesh/CurvilinearMesh.py
similarity index 97%
rename from SimPEG/Mesh/LogicallyRectMesh.py
rename to SimPEG/Mesh/CurvilinearMesh.py
index b2dc6f55..f8b0bcd2 100644
--- a/SimPEG/Mesh/LogicallyRectMesh.py
+++ b/SimPEG/Mesh/CurvilinearMesh.py
@@ -10,24 +10,24 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
 normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
 
 
-class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
+class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
     """
-    LogicallyRectMesh is a mesh class that deals with logically rectangular meshes.
+    CurvilinearMesh is a mesh class that deals with curvilinear meshes.
 
-    Example of a logically rectangular mesh:
+    Example of a curvilinear mesh:
 
     .. plot::
             :include-source:
 
             from SimPEG import Mesh, Utils
             X, Y = Utils.exampleLrmGrid([3,3],'rotate')
-            M = Mesh.LogicallyRectMesh([X, Y])
+            M = Mesh.CurvilinearMesh([X, Y])
             M.plotGrid(showIt=True)
     """
 
     __metaclass__ = Utils.SimPEGMetaClass
 
-    _meshType = 'LRM'
+    _meshType = 'Curv'
 
     def __init__(self, nodes):
         assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray"
@@ -38,7 +38,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
             assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % i)
 
         assert len(nodes[0].shape) == len(nodes), "Dimension mismatch"
-        assert len(nodes[0].shape) > 1, "Not worth using LRM for a 1D mesh."
+        assert len(nodes[0].shape) > 1, "Not worth using Curv for a 1D mesh."
 
         BaseRectangularMesh.__init__(self, np.array(nodes[0].shape)-1, None)
 
@@ -343,7 +343,7 @@ class LogicallyRectMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
 
             from SimPEG import Mesh, Utils
             X, Y = Utils.exampleLrmGrid([3,3],'rotate')
-            M = Mesh.LogicallyRectMesh([X, Y])
+            M = Mesh.CurvilinearMesh([X, Y])
             M.plotGrid(showIt=True)
 
         """
@@ -435,9 +435,9 @@ if __name__ == '__main__':
     dee3 = True
     if dee3:
         X, Y, Z = Utils.ndgrid(h1, h2, h3, vector=False)
-        M = LogicallyRectMesh([X, Y, Z])
+        M = CurvilinearMesh([X, Y, Z])
     else:
         X, Y = Utils.ndgrid(h1, h2, vector=False)
-        M = LogicallyRectMesh([X, Y])
+        M = CurvilinearMesh([X, Y])
 
     print M.r(M.normals, 'F', 'Fx', 'V')
diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
index 2cf48a7c..42a08702 100644
--- a/SimPEG/Mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -328,7 +328,7 @@ class InnerProducts(object):
         iijj = ndgrid(i, j)
         ii, jj = iijj[:, 0], iijj[:, 1]
 
-        if M._meshType == 'LRM':
+        if M._meshType == 'Curv':
             fN1 = M.r(M.normals, 'F', 'Fx', 'M')
             fN2 = M.r(M.normals, 'F', 'Fy', 'M')
 
@@ -353,7 +353,7 @@ class InnerProducts(object):
 
             PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF))
 
-            if M._meshType == 'LRM':
+            if M._meshType == 'Curv':
                 I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]),
                                            getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy]))
                 PXX = I2x2 * PXX
@@ -376,7 +376,7 @@ class InnerProducts(object):
         iijjkk = ndgrid(i, j, k)
         ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
 
-        if M._meshType == 'LRM':
+        if M._meshType == 'Curv':
             fN1 = M.r(M.normals, 'F', 'Fx', 'M')
             fN2 = M.r(M.normals, 'F', 'Fy', 'M')
             fN3 = M.r(M.normals, 'F', 'Fz', 'M')
@@ -410,7 +410,7 @@ class InnerProducts(object):
 
             PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr()
 
-            if M._meshType == 'LRM':
+            if M._meshType == 'Curv':
                 I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]),
                                            getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]),
                                            getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ]))
@@ -432,7 +432,7 @@ class InnerProducts(object):
         iijj = ndgrid(i, j)
         ii, jj = iijj[:, 0], iijj[:, 1]
 
-        if M._meshType == 'LRM':
+        if M._meshType == 'Curv':
             eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
             eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
 
@@ -452,7 +452,7 @@ class InnerProducts(object):
 
             PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr()
 
-            if M._meshType == 'LRM':
+            if M._meshType == 'Curv':
                 I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]),
                                            getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j]))
                 PXX = I2x2 * PXX
@@ -466,7 +466,7 @@ class InnerProducts(object):
         iijjkk = ndgrid(i, j, k)
         ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
 
-        if M._meshType == 'LRM':
+        if M._meshType == 'Curv':
             eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
             eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
             eT3 = M.r(M.tangents, 'E', 'Ez', 'M')
@@ -495,7 +495,7 @@ class InnerProducts(object):
 
             PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr()
 
-            if M._meshType == 'LRM':
+            if M._meshType == 'Curv':
                 I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]),
                                            getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]),
                                            getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k]))
diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py
index b5ede21f..38432d3c 100644
--- a/SimPEG/Mesh/View.py
+++ b/SimPEG/Mesh/View.py
@@ -328,6 +328,7 @@ class TensorView(object):
             v = getattr(np,view)(v) # e.g. np.real(v)
             if clim is None:
                 clim = [v.min(),v.max()]
+            v = np.ma.masked_where(np.isnan(v), v)
             out += (ax.pcolormesh(self.vectorNx, self.vectorNy, v.T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
         elif view in ['vec']:
             U, V = self.r(v.reshape((self.nC,-1), order='F'), 'CC', 'CC', 'M')
@@ -576,11 +577,11 @@ class CylView(object):
     def plotImage(self, *args, **kwargs):
         return self._plotCylTensorMesh('plotImage', *args, **kwargs)
 
-class LomView(object):
+class CurvView(object):
     """
-    Provides viewing functions for LogicallyOrthogonalMesh
+    Provides viewing functions for CurvilinearMesh
 
-    This class is inherited by LogicallyOrthogonalMesh
+    This class is inherited by CurvilinearMesh
 
     """
     def __init__(self):
@@ -594,8 +595,8 @@ class LomView(object):
             :include-source:
 
             from SimPEG import Mesh, Utils
-            X, Y = Utils.exampleLomGird([3,3],'rotate')
-            M = Mesh.LogicallyOrthogonalMesh([X, Y])
+            X, Y = Utils.exampleCurvGird([3,3],'rotate')
+            M = Mesh.CurvilinearMesh([X, Y])
             M.plotGrid(showIt=True)
 
         """
diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py
index deb11321..9aadfcba 100644
--- a/SimPEG/Mesh/__init__.py
+++ b/SimPEG/Mesh/__init__.py
@@ -1,5 +1,5 @@
 from TensorMesh import TensorMesh
 from CylMesh import CylMesh
-from LogicallyRectMesh import LogicallyRectMesh
+from CurvilinearMesh import CurvilinearMesh
 from TreeMesh import TreeMesh
 from BaseMesh import BaseMesh
diff --git a/SimPEG/Models.py b/SimPEG/Models.py
index f69b749c..65cfdeee 100644
--- a/SimPEG/Models.py
+++ b/SimPEG/Models.py
@@ -31,5 +31,5 @@ class Model(np.ndarray):
     @property
     def transformDeriv(self):
         if getattr(self, '_transformDeriv', None) is None:
-            self.deriv = self.mapping.deriv(self.view(np.ndarray))
-        return self.deriv
+            self._transformDeriv = self.mapping.deriv(self.view(np.ndarray))
+        return self._transformDeriv
diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py
index 549cb8ae..4f2cb062 100644
--- a/SimPEG/Optimization.py
+++ b/SimPEG/Optimization.py
@@ -97,6 +97,8 @@ class Minimize(object):
     tolG = 1e-1        #: Tolerance on gradient norm
     eps = 1e-5         #: Small value
 
+    stopNextIteration = False #: Stops the optimization program nicely.
+
     debug   = False    #: Print debugging information
     debugLS = False    #: Print debugging information for the line-search
 
@@ -186,6 +188,7 @@ class Minimize(object):
                 xt, caught = self.modifySearchDirectionBreak(p)
                 if not caught: return self.xc
             self.doEndIteration(xt)
+            if self.stopNextIteration: break
 
         self.printDone()
         self.finish()
@@ -210,6 +213,7 @@ class Minimize(object):
 
         self.iter = 0
         self.iterLS = 0
+        self.stopNextIteration = False
 
         x0 = self.projection(x0)  # ensure that we start of feasible.
         self.x0 = x0
diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py
index dff99a97..607cff5b 100644
--- a/SimPEG/Problem.py
+++ b/SimPEG/Problem.py
@@ -1,279 +1,7 @@
 import Utils, Survey, Models, numpy as np, scipy.sparse as sp
 Solver = Utils.SolverUtils.Solver
 import Maps, Mesh
-
-
-class Fields(object):
-    """Fancy Field Storage
-
-        u[:,'phi'] = phi
-        print u[src0,'phi']
-
-    """
-
-    knownFields = None  #: Known fields,   a dict with locations,         e.g. {"e": "E", "phi": "CC"}
-    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.
-
-    def __init__(self, mesh, survey, **kwargs):
-        self.survey = survey
-        self.mesh = mesh
-        Utils.setKwargs(self, **kwargs)
-        self._fields = {}
-
-        if self.knownFields is None:
-            raise Exception('knownFields cannot be set to None')
-        if self.aliasFields is None:
-            self.aliasFields = {}
-
-        allFields = [k for k in self.knownFields] + [a for a in self.aliasFields]
-        assert len(allFields) == len(set(allFields)), 'Aliased fields and Known Fields have overlapping definitions.'
-        self.startup()
-
-    def startup(self):
-        pass
-
-    @property
-    def approxSize(self):
-        """The approximate cost to storing all of the known fields."""
-        sz = 0.0
-        for f in self.knownFields:
-            loc =self.knownFields[f]
-            sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2)
-        return "%e MB"%sz
-
-    def _storageShape(self, loc):
-        nSrc = self.survey.nSrc
-
-        nP = {'CC': self.mesh.nC,
-              'N':  self.mesh.nN,
-              'F':  self.mesh.nF,
-              'E':  self.mesh.nE}[loc]
-
-        return (nP, nSrc)
-
-    def _initStore(self, name):
-        if name in self._fields:
-            return self._fields[name]
-
-        assert name in self.knownFields, 'field name is not known.'
-
-        loc = self.knownFields[name]
-
-        if type(self.dtype) is dict:
-            dtype = self.dtype[name]
-        else:
-            dtype = self.dtype
-        field = np.zeros(self._storageShape(loc), dtype=dtype)
-
-        self._fields[name] = field
-
-        return field
-
-    def _srcIndex(self, srcTestList):
-        if type(srcTestList) is slice:
-            ind = srcTestList
-        else:
-            if type(srcTestList) is not list:
-                srcTestList = [srcTestList]
-            for srcTest in srcTestList:
-                if srcTest not in self.survey.srcList:
-                    raise KeyError('Invalid Source, not in survey list.')
-
-            ind = np.in1d(self.survey.srcList, srcTestList)
-        return ind
-
-    def _nameIndex(self, name, accessType):
-
-        if type(name) is slice:
-            assert name == slice(None,None,None), 'Fancy field name slicing is not supported... yet.'
-            name = None
-
-        if name is None:
-            return
-        if accessType=='set' and name not in self.knownFields:
-            if name in self.aliasFields:
-                raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name)
-            else:
-                raise KeyError('Invalid field name (%s) for setter'%name)
-
-        elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields):
-            raise KeyError('Invalid field name (%s) for getter'%name)
-        return name
-
-    def _indexAndNameFromKey(self, key, accessType):
-        if type(key) is not tuple:
-            key = (key,)
-        if len(key) == 1:
-            key += (None,)
-
-        assert len(key) == 2, 'must be [Src, fieldName]'
-
-        srcTestList, name = key
-        name = self._nameIndex(name, accessType)
-        ind = self._srcIndex(srcTestList)
-        return ind, name
-
-    def __setitem__(self, key, value):
-        ind, name = self._indexAndNameFromKey(key, 'set')
-        if name is None:
-            freq = key
-            assert type(value) is dict, 'New fields must be a dictionary, if field is not specified.'
-            newFields = value
-        elif name in self.knownFields:
-            newFields = {name: value}
-        else:
-            raise Exception('Unknown setter')
-
-        for name in newFields:
-            field = self._initStore(name)
-            self._setField(field, newFields[name], name, ind)
-
-    def __getitem__(self, key):
-        ind, name = self._indexAndNameFromKey(key, 'get')
-        if name is None:
-            out = {}
-            for name in self._fields:
-                out[name] = self._getField(name, ind)
-            return out
-        return self._getField(name, ind)
-
-    def _setField(self, field, val, name, ind):
-        if isinstance(val, np.ndarray) and (field.shape[1] == 1 or val.ndim == 1):
-            val = Utils.mkvc(val,2)
-        field[:,ind] = val
-
-    def _getField(self, name, ind):
-        if name in self._fields:
-            out = self._fields[name][:,ind]
-        else:
-            # Aliased fields
-            alias, loc, func = self.aliasFields[name]
-
-            srcII   = np.array(self.survey.srcList)[ind]
-            if isinstance(srcII, np.ndarray):
-                srcII = srcII.tolist()
-            if len(srcII) == 1:
-                srcII = srcII[0]
-
-            if type(func) is str:
-                assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
-                func = getattr(self, func)
-            out = func(self._fields[alias][:,ind], srcII)
-
-        if out.shape[1] == 1:
-            out = Utils.mkvc(out)
-        return out
-
-    def __contains__(self, other):
-        if other in self.aliasFields:
-            other = self.aliasFields[other][0]
-        return self._fields.__contains__(other)
-
-
-class TimeFields(Fields):
-    """Fancy Field Storage for time domain problems
-
-        u[:,'phi', timeInd] = phi
-        print u[src0,'phi']
-
-    """
-
-    def _storageShape(self, loc):
-        nP = {'CC': self.mesh.nC,
-              'N':  self.mesh.nN,
-              'F':  self.mesh.nF,
-              'E':  self.mesh.nE}[loc]
-        nSrc = self.survey.nSrc
-        nT = self.survey.prob.nT + 1
-        return (nP, nSrc, nT)
-
-    def _indexAndNameFromKey(self, key, accessType):
-        if type(key) is not tuple:
-            key = (key,)
-        if len(key) == 1:
-            key += (None,)
-        if len(key) == 2:
-            key += (slice(None,None,None),)
-
-        assert len(key) == 3, 'must be [Src, fieldName, times]'
-
-        srcTestList, name, timeInd = key
-
-        name = self._nameIndex(name, accessType)
-        srcInd = self._srcIndex(srcTestList)
-
-        return (srcInd, timeInd), name
-
-    def _correctShape(self, name, ind, deflate=False):
-        srcInd, timeInd = ind
-        if name in self.knownFields:
-            loc = self.knownFields[name]
-        else:
-            loc = self.aliasFields[name][1]
-        nP, total_nSrc, total_nT = self._storageShape(loc)
-        nSrc = np.ones(total_nSrc, dtype=bool)[srcInd].sum()
-        nT  = np.ones(total_nT, dtype=bool)[timeInd].sum()
-        shape = nP, nSrc, nT
-        if deflate:
-             shape = tuple([s for s in shape if s > 1])
-        return shape
-
-    def _setField(self, field, val, name, ind):
-        srcInd, timeInd = ind
-        shape = self._correctShape(name, ind)
-        if Utils.isScalar(val):
-            field[:,srcInd,timeInd] = val
-            return
-        if val.size != np.array(shape).prod():
-            raise ValueError('Incorrect size for data.')
-        correctShape = field[:,srcInd,timeInd].shape
-        field[:,srcInd,timeInd] = val.reshape(correctShape, order='F')
-
-    def _getField(self, name, ind):
-        srcInd, timeInd = ind
-
-        if name in self._fields:
-            out = self._fields[name][:,srcInd,timeInd]
-        else:
-            # Aliased fields
-            alias, loc, func = self.aliasFields[name]
-            if type(func) is str:
-                assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
-                func = getattr(self, func)
-            pointerFields = self._fields[alias][:,srcInd,timeInd]
-            pointerShape = self._correctShape(alias, ind)
-            pointerFields = pointerFields.reshape(pointerShape, order='F')
-
-            timeII = np.arange(self.survey.prob.nT + 1)[timeInd]
-            srcII   = np.array(self.survey.srcList)[srcInd]
-            if isinstance(srcII, np.ndarray):
-                srcII = srcII.tolist()
-            if len(srcII) == 1:
-                srcII = srcII[0]
-
-            if timeII.size == 1:
-                pointerShapeDeflated = self._correctShape(alias, ind, deflate=True)
-                pointerFields = pointerFields.reshape(pointerShapeDeflated, order='F')
-                out = func(pointerFields, srcII, timeII)
-            else: #loop over the time steps
-                nT = pointerShape[2]
-                out = range(nT)
-                for i, TIND_i in enumerate(timeII):
-                    fieldI = pointerFields[:,:,i]
-                    if fieldI.ndim == 2 and fieldI.shape[1] == 1:
-                        fieldI = Utils.mkvc(fieldI)
-                    out[i] = func(fieldI, srcII, TIND_i)
-                    if out[i].ndim == 1:
-                        out[i] = out[i][:,np.newaxis,np.newaxis]
-                    elif out[i].ndim == 2:
-                        out[i] = out[i][:,:,np.newaxis]
-                out = np.concatenate(out, axis=2)
-
-        shape = self._correctShape(name, ind, deflate=True)
-        return out.reshape(shape, order='F')
-
-
+from Fields import Fields, TimeFields
 
 class BaseProblem(object):
     """
@@ -290,15 +18,27 @@ class BaseProblem(object):
     Solver = Solver   #: A SimPEG Solver class.
     solverOpts = {}   #: Sovler options as a kwarg dict
 
-    mapping = None    #: A SimPEG.Map instance.
     mesh    = None    #: A SimPEG.Mesh instance.
 
+    PropMap = None    #: A SimPEG PropertyMap class.
+
+    @property
+    def mapping(self):
+        "A SimPEG.Map instance or a property map is PropMap is not None"
+        return getattr(self, '_mapping', None)
+    @mapping.setter
+    def mapping(self, val):
+        if self.PropMap is None:
+            val._assertMatchesPair(self.mapPair)
+            self._mapping = val
+        else:
+            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."
         self.mesh = mesh
         self.mapping = mapping or Maps.IdentityMap(mesh)
-        self.mapping._assertMatchesPair(self.mapPair)
 
     @property
     def survey(self):
@@ -334,7 +74,10 @@ class BaseProblem(object):
     def curModel(self, value):
         if value is self.curModel:
             return # it is the same!
-        self._curModel = Models.Model(value, self.mapping)
+        if self.PropMap is not None:
+            self._curModel = self.mapping(value)
+        else:
+            self._curModel = Models.Model(value, self.mapping)
         for prop in self.deleteTheseOnModelUpdate:
             if hasattr(self, prop):
                 delattr(self, prop)
diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py
new file mode 100644
index 00000000..995216f7
--- /dev/null
+++ b/SimPEG/PropMaps.py
@@ -0,0 +1,260 @@
+import Utils, Maps, numpy as np, scipy.sparse as sp
+
+class Property(object):
+
+    name           = ''
+    doc            = ''
+
+    defaultVal     = None
+    defaultInvProp = False
+
+    def __init__(self, doc, **kwargs):
+        # Set the default after all other params are set
+        self.doc = doc
+        Utils.setKwargs(self, **kwargs)
+
+    @property
+    def propertyLink(self):
+        "Can be something like: ('sigma', Maps.ReciprocalMap)"
+        return getattr(self, '_propertyLink', None)
+    @propertyLink.setter
+    def propertyLink(self, value):
+        assert type(value) is tuple and len(value) == 2 and type(value[0]) is str and issubclass(value[1], Maps.IdentityMap), 'Use format: ("%s", Maps.ReciprocalMap)'%self.name
+        self._propertyLink = value
+
+    def _getMapProperty(self):
+        prop = self
+        def fget(self):
+            return getattr(self, '_%sMap'%prop.name, None)
+        def fset(self, val):
+            if prop.propertyLink is not None:
+                linkName, linkMap = prop.propertyLink
+                assert getattr(self, '%sMap'%linkName, None) is None, 'Cannot set both sides of a linked property.'
+            # TODO: Check if the mapping can be correct
+            setattr(self, '_%sMap'%prop.name, val)
+        return property(fget=fget, fset=fset, doc=prop.doc)
+
+    def _getIndexProperty(self):
+        prop = self
+        def fget(self):
+            return getattr(self, '_%sIndex'%prop.name, slice(None))
+        def fset(self, val):
+            setattr(self, '_%sIndex'%prop.name, val)
+        return property(fget=fget, fset=fset, doc=prop.doc)
+
+    def _getProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None and prop.propertyLink is None:
+                return prop.defaultVal
+
+            if mapping is None and prop.propertyLink is not None:
+                linkName, linkMapClass = prop.propertyLink
+                linkMap = linkMapClass(None)
+                if getattr(self, '%sMap'%linkName, None) is None:
+                    return prop.defaultVal
+                m = getattr(self, '%s'%linkName)
+                return linkMap * m
+
+            m = getattr(self, '%sModel'%prop.name)
+            return mapping * m
+        return property(fget=fget)
+
+    def _getModelDerivProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None and prop.propertyLink is None:
+                return None
+
+            if mapping is None and prop.propertyLink is not None:
+                linkName, linkMapClass = prop.propertyLink
+                linkedMap = getattr(self, '%sMap'%linkName)
+                if linkedMap is None:
+                    return None
+                linkMap = linkMapClass(None) * linkedMap
+                m = getattr(self, '%s'%linkName)
+                return linkMap.deriv( m )
+
+            m = getattr(self, '%sModel'%prop.name)
+            return mapping.deriv( m )
+        return property(fget=fget)
+
+    def _getModelProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None:
+                return None
+            index = getattr(self.propMap, '%sIndex'%prop.name)
+            return self.vector[index]
+        return property(fget=fget)
+
+    def _getModelProjProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None:
+                return None
+            inds = getattr(self.propMap, '%sIndex'%prop.name)
+            if type(inds) is slice:
+                inds = range(*inds.indices(self.nP))
+            nI, nP = len(inds),self.nP
+            return sp.csr_matrix((np.ones(nI), (range(nI), inds) ), shape=(nI, nP))
+        return property(fget=fget)
+
+    def _getModelMapProperty(self):
+        prop = self
+        def fget(self):
+            return getattr(self.propMap, '_%sMap'%prop.name, None)
+        return property(fget=fget)
+
+
+
+class PropModel(object):
+    def __init__(self, propMap, vector):
+        self.propMap = propMap
+        self.vector  = vector
+        assert len(self.vector) == self.nP
+
+    @property
+    def nP(self):
+        inds = []
+        if getattr(self, '_nP', None) is None:
+            for name in self.propMap._properties:
+                index = getattr(self.propMap, '%sIndex'%name, None)
+                if index is not None:
+                    if type(index) is slice:
+                        inds += range(*index.indices(len(self.vector)))
+                    else:
+                        inds += list(index)
+            self._nP = len(set(inds))
+        return self._nP
+
+    def __contains__(self, val):
+        return val in self.propMap
+
+
+
+_PROPMAPCLASSREGISTRY = {}
+
+class _PropMapMetaClass(type):
+    def __new__(cls, name, bases, attrs):
+        assert name.endswith('PropMap'), 'Please use convention: ___PropMap, e.g. ElectromagneticPropMap'
+        _properties = {}
+        for base in bases:
+            for baseProp in getattr(base, '_properties', {}):
+                _properties[baseProp] = base._properties[baseProp]
+        keys = [key for key in attrs]
+        for attr in keys:
+            if isinstance(attrs[attr], Property):
+                attrs[attr].name = attr
+                attrs[attr + 'Map'  ] = attrs[attr]._getMapProperty()
+                attrs[attr + 'Index'] = attrs[attr]._getIndexProperty()
+                _properties[attr] = attrs[attr]
+                attrs.pop(attr)
+
+        attrs['_properties'] = _properties
+
+        defaultInvProps = []
+        for p in _properties:
+            prop = _properties[p]
+            if prop.defaultInvProp:
+                defaultInvProps += [p]
+            if prop.propertyLink is not None:
+                assert prop.propertyLink[0] in _properties, "You can only link to things that exist: '%s' is trying to link to '%s'"%(prop.name, prop.propertyLink[0])
+        if len(defaultInvProps) > 1:
+            raise Exception('You have more than one default inversion property: %s' % defaultInvProps)
+
+        newClass = super(_PropMapMetaClass, cls).__new__(cls, name, bases, attrs)
+
+        newClass.PropModel = cls.createPropModelClass(newClass, name, _properties)
+
+        _PROPMAPCLASSREGISTRY[name] = newClass
+        return newClass
+
+    def createPropModelClass(self, name, _properties):
+
+        attrs = dict()
+
+        for attr in _properties:
+            prop = _properties[attr]
+
+            attrs[attr          ] = prop._getProperty()
+            attrs[attr + 'Map'  ] = prop._getModelMapProperty()
+            attrs[attr + 'Proj' ] = prop._getModelProjProperty()
+            attrs[attr + 'Model'] = prop._getModelProperty()
+            attrs[attr + 'Deriv'] = prop._getModelDerivProperty()
+
+        return type(name.replace('PropMap', 'PropModel'), (PropModel, ), attrs)
+
+
+class PropMap(object):
+    __metaclass__ = _PropMapMetaClass
+
+    def __init__(self, mappings):
+        """
+            PropMap takes a multi parameter model and maps it to the equivalent PropModel
+        """
+        if type(mappings) is dict:
+            assert np.all([k in ['maps', 'slices'] for k in mappings]), 'Dict must only have properties "maps" and "slices"'
+            self.setup(mappings['maps'], slices=mappings['slices'])
+        elif type(mappings) is list:
+            self.setup(mappings)
+        elif isinstance(mappings, Maps.IdentityMap):
+            self.setup([(self.defaultInvProp, mappings)])
+        else:
+            raise Exception('mappings must be a dict, a mapping, or a list of tuples.')
+
+
+    def setup(self, maps, slices=None):
+        """
+            Sets up the maps and slices for the PropertyMap
+
+
+            :param list maps: [('sigma', sigmaMap), ('mu', muMap), ...]
+            :param list slices: [('sigma', slice(0,nP)), ('mu', [1,2,5,6]), ...]
+
+        """
+        assert np.all([
+                type(m) is tuple and
+                len(m)==2 and
+                type(m[0]) is str and
+                m[0] in self._properties and
+                isinstance(m[1], Maps.IdentityMap)
+                for m in maps]), "Use signature: [%s]" % (', '.join(["('%s', %sMap)"%(p,p) for p in self._properties]))
+        if slices is None:
+            slices = dict()
+        else:
+            assert np.all([
+                s in self._properties and
+                (type(slices[s]) in [slice, list] or isinstance(slices[s], np.ndarray))
+                for s in slices]), 'Slices must be for each property'
+
+        self.clearMaps()
+
+        nP = 0
+        for name, mapping in maps:
+            setattr(self, '%sMap'%name, mapping)
+            setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP)))
+            nP += mapping.nP
+
+    @property
+    def defaultInvProp(self):
+        for name in self._properties:
+            p = self._properties[name]
+            if p.defaultInvProp:
+                return p.name
+
+    def clearMaps(self):
+        for name in self._properties:
+            setattr(self, '%sMap'%name, None)
+            setattr(self, '%sIndex'%name, None)
+
+    def __call__(self, vec):
+        return self.PropModel(self, vec)
+
+    def __contains__(self, val):
+        activeMaps = [name for name in self._properties if getattr(self, '%sMap'%name) is not None]
+        return val in activeMaps
diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py
index 10e4c76f..0fdb0cd1 100644
--- a/SimPEG/Survey.py
+++ b/SimPEG/Survey.py
@@ -1,4 +1,4 @@
-import Utils, numpy as np, scipy.sparse as sp
+import Utils, numpy as np, scipy.sparse as sp, uuid
 
 
 class BaseRx(object):
@@ -13,6 +13,7 @@ class BaseRx(object):
     storeProjections = True #: Store calls to getP (organized by mesh)
 
     def __init__(self, locs, rxType, **kwargs):
+        self.uid = str(uuid.uuid4())
         self.locs = locs
         self.rxType = rxType
         self._Ps = {}
@@ -119,14 +120,12 @@ class BaseSrc(object):
     rxList = None #: SimPEG Receiver List
     rxPair = BaseRx
 
-    def __init__(self, loc, srcType, rxList, **kwargs):
+    def __init__(self, rxList, **kwargs):
         assert type(rxList) is list, 'rxList must be a list'
         for rx in rxList:
             assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxPair.__name__
         assert len(set(rxList)) == len(rxList), 'The rxList must be unique'
-
-        self.loc    = loc
-        self.srcType = srcType
+        self.uid = str(uuid.uuid4())
         self.rxList = rxList
         Utils.setKwargs(self, **kwargs)
 
@@ -146,6 +145,7 @@ class Data(object):
     """Fancy data storage by Src and Rx"""
 
     def __init__(self, survey, v=None):
+        self.uid = str(uuid.uuid4())
         self.survey = survey
         self._dataDict = {}
         for src in self.survey.srcList:
@@ -227,6 +227,19 @@ class BaseSurvey(object):
         assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of %s' % self.srcPair.__name__
         assert len(set(value)) == len(value), 'The srcList must be unique'
         self._srcList = value
+        self._sourceOrder = dict()
+        [self._sourceOrder.setdefault(src.uid, ii) for ii, src in enumerate(self._srcList)]
+
+    def getSourceIndex(self, sources):
+        if type(sources) is not list:
+            sources = [sources]
+        for src in sources:
+            if getattr(src,'uid',None) is None:
+                raise KeyError('Source does not have a uid: %s'%str(src))
+        inds = map(lambda src: self._sourceOrder.get(src.uid, None), sources)
+        if None in inds:
+            raise KeyError('Some of the sources specified are not in this survey. %s'%str(inds))
+        return inds
 
     @property
     def prob(self):
@@ -360,3 +373,4 @@ class BaseSurvey(object):
         noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape)
         self.dobs = self.dtrue+noise
         self.std = self.dobs*0 + std
+        return self.dobs
diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py
index dbcb92b6..c6aa5bc4 100644
--- a/SimPEG/Tests/TestUtils.py
+++ b/SimPEG/Tests/TestUtils.py
@@ -3,7 +3,7 @@ import matplotlib.pyplot as plt
 from numpy.linalg import norm
 from SimPEG.Utils import mkvc, sdiag, diagEst
 from SimPEG import Utils
-from SimPEG.Mesh import TensorMesh, LogicallyRectMesh, CylMesh
+from SimPEG.Mesh import TensorMesh, CurvilinearMesh, CylMesh
 import numpy as np
 import scipy.sparse as sp
 import unittest
@@ -115,7 +115,7 @@ class OrderTest(unittest.TestCase):
                 max_h = max([np.max(hi) for hi in self.M.h])
             return max_h
 
-        elif 'LRM' in self._meshType:
+        elif 'Curv' in self._meshType:
             if 'uniform' in self._meshType:
                 kwrd = 'rect'
             elif 'rotate' in self._meshType:
@@ -126,10 +126,10 @@ class OrderTest(unittest.TestCase):
                 raise Exception('Lom not supported for 1D')
             elif self.meshDimension == 2:
                 X, Y = Utils.exampleLrmGrid([nc, nc], kwrd)
-                self.M = LogicallyRectMesh([X, Y])
+                self.M = CurvilinearMesh([X, Y])
             elif self.meshDimension == 3:
                 X, Y, Z = Utils.exampleLrmGrid([nc, nc, nc], kwrd)
-                self.M = LogicallyRectMesh([X, Y, Z])
+                self.M = CurvilinearMesh([X, Y, Z])
             return 1./nc
 
     def getError(self):
diff --git a/SimPEG/Tests/test_CurvilinearMesh.py b/SimPEG/Tests/test_CurvilinearMesh.py
new file mode 100644
index 00000000..42e3d877
--- /dev/null
+++ b/SimPEG/Tests/test_CurvilinearMesh.py
@@ -0,0 +1,104 @@
+import numpy as np
+import unittest
+from SimPEG.Mesh import TensorMesh, CurvilinearMesh
+from SimPEG.Utils import ndgrid
+
+
+class BasicCurvTests(unittest.TestCase):
+
+    def setUp(self):
+        a = np.array([1, 1, 1])
+        b = np.array([1, 2])
+        c = np.array([1, 4])
+        gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h]
+        X, Y = ndgrid(gridIt([a, b]), vector=False)
+        self.TM2 = TensorMesh([a, b])
+        self.Curv2 = CurvilinearMesh([X, Y])
+        X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False)
+        self.TM3 = TensorMesh([a, b, c])
+        self.Curv3 = CurvilinearMesh([X, Y, Z])
+
+    def test_area_3D(self):
+        test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2])
+        self.assertTrue(np.all(self.Curv3.area == test_area))
+
+    def test_vol_3D(self):
+        test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8])
+        np.testing.assert_almost_equal(self.Curv3.vol, test_vol)
+        self.assertTrue(True)  # Pass if you get past the assertion.
+
+    def test_vol_2D(self):
+        test_vol = np.array([1, 1, 1, 2, 2, 2])
+        t1 = np.all(self.Curv2.vol == test_vol)
+        self.assertTrue(t1)
+
+    def test_edge_3D(self):
+        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
+        t1 = np.all(self.Curv3.edge == test_edge)
+        self.assertTrue(t1)
+
+    def test_edge_2D(self):
+        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2])
+        t1 = np.all(self.Curv2.edge == test_edge)
+        self.assertTrue(t1)
+
+    def test_tangents(self):
+        T = self.Curv2.tangents
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv2.nEx)))
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv2.nEx)))
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv2.nEy)))
+        self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv2.nEy)))
+
+        T = self.Curv3.tangents
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv3.nEx)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv3.nEx)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.Curv3.nEx)))
+
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv3.nEy)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv3.nEy)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.Curv3.nEy)))
+
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.Curv3.nEz)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.Curv3.nEz)))
+        self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.Curv3.nEz)))
+
+    def test_normals(self):
+        N = self.Curv2.normals
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv2.nFx)))
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv2.nFx)))
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv2.nFy)))
+        self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv2.nFy)))
+
+        N = self.Curv3.normals
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv3.nFx)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv3.nFx)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.Curv3.nFx)))
+
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv3.nFy)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv3.nFy)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.Curv3.nFy)))
+
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.Curv3.nFz)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.Curv3.nFz)))
+        self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.Curv3.nFz)))
+
+    def test_grid(self):
+        self.assertTrue(np.all(self.Curv2.gridCC == self.TM2.gridCC))
+        self.assertTrue(np.all(self.Curv2.gridN == self.TM2.gridN))
+        self.assertTrue(np.all(self.Curv2.gridFx == self.TM2.gridFx))
+        self.assertTrue(np.all(self.Curv2.gridFy == self.TM2.gridFy))
+        self.assertTrue(np.all(self.Curv2.gridEx == self.TM2.gridEx))
+        self.assertTrue(np.all(self.Curv2.gridEy == self.TM2.gridEy))
+
+        self.assertTrue(np.all(self.Curv3.gridCC == self.TM3.gridCC))
+        self.assertTrue(np.all(self.Curv3.gridN == self.TM3.gridN))
+        self.assertTrue(np.all(self.Curv3.gridFx == self.TM3.gridFx))
+        self.assertTrue(np.all(self.Curv3.gridFy == self.TM3.gridFy))
+        self.assertTrue(np.all(self.Curv3.gridFz == self.TM3.gridFz))
+        self.assertTrue(np.all(self.Curv3.gridEx == self.TM3.gridEx))
+        self.assertTrue(np.all(self.Curv3.gridEy == self.TM3.gridEy))
+        self.assertTrue(np.all(self.Curv3.gridEz == self.TM3.gridEz))
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/SimPEG/Tests/test_Survey.py b/SimPEG/Tests/test_Fields.py
similarity index 81%
rename from SimPEG/Tests/test_Survey.py
rename to SimPEG/Tests/test_Fields.py
index f93395ce..629189b8 100644
--- a/SimPEG/Tests/test_Survey.py
+++ b/SimPEG/Tests/test_Fields.py
@@ -1,7 +1,8 @@
 import unittest
 from SimPEG import *
 
-class DataAndFieldsTest(unittest.TestCase):
+
+class FieldsTest(unittest.TestCase):
 
     def setUp(self):
         mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
@@ -9,14 +10,14 @@ class DataAndFieldsTest(unittest.TestCase):
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
         srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
-        Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
+        Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
+        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
         srcList = [Src0,Src1,Src2,Src3,Src4]
         survey = Survey.BaseSurvey(srcList=srcList)
         self.D = Survey.Data(survey)
@@ -26,25 +27,6 @@ class DataAndFieldsTest(unittest.TestCase):
         self.mesh = mesh
         self.XYZ = XYZ
 
-    def test_overlappingFields(self):
-        self.assertRaises(AssertionError, Problem.Fields, self.F.mesh, self.F.survey,
-                            knownFields={'b':'F'},
-                            aliasFields={'b':['b',(lambda F, b, ind: b)]})
-
-    def test_data(self):
-        V = []
-        for src in self.D.survey.srcList:
-            for rx in src.rxList:
-                v = np.random.rand(rx.nD)
-                V += [v]
-                self.D[src, rx] = v
-                self.assertTrue(np.all(v == self.D[src, rx]))
-        V = np.concatenate(V)
-        self.assertTrue(np.all(V == Utils.mkvc(self.D)))
-
-        D2 = Survey.Data(self.D.survey, V)
-        self.assertTrue(np.all(Utils.mkvc(D2) == Utils.mkvc(self.D)))
-
     def test_contains(self):
         F = self.F
         nSrc = F.survey.nSrc
@@ -55,10 +37,10 @@ class DataAndFieldsTest(unittest.TestCase):
         self.assertTrue('b' not in F)
         self.assertTrue('e' in F)
 
-    def test_uniqueSrcs(self):
-        srcs = self.D.survey.srcList
-        srcs += [srcs[0]]
-        self.assertRaises(AssertionError, Survey.BaseSurvey, srcList=srcs)
+    def test_overlappingFields(self):
+        self.assertRaises(AssertionError, Problem.Fields, self.F.mesh, self.F.survey,
+                            knownFields={'b':'F'},
+                            aliasFields={'b':['b',(lambda F, b, ind: b)]})
 
     def test_SetGet(self):
         F = self.F
@@ -80,9 +62,9 @@ class DataAndFieldsTest(unittest.TestCase):
 
         b = np.random.rand(F.mesh.nF,1)
         F[self.Src0, 'b'] = b
-        self.assertTrue(np.all(F[self.Src0, 'b'] == Utils.mkvc(b)))
+        self.assertTrue(np.all(F[self.Src0, 'b'] == b))
 
-        b = np.random.rand(F.mesh.nF)
+        b = np.random.rand(F.mesh.nF,1)
         F[self.Src0, 'b'] = b
         self.assertTrue(np.all(F[self.Src0, 'b'] == b))
 
@@ -96,10 +78,10 @@ class DataAndFieldsTest(unittest.TestCase):
 
         b = np.random.rand(F.mesh.nF, 2)
         F[[self.Src0, self.Src1],'b'] = b
-        self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,))
-        self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,))
-        self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0]))
-        self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1]))
+        self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,1))
+        self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,1))
+        self.assertTrue(np.all(F[self.Src0,'b'] == Utils.mkvc(b[:,0],2)))
+        self.assertTrue(np.all(F[self.Src1,'b'] == Utils.mkvc(b[:,1],2)))
 
     def test_assertions(self):
         freq = [self.Src0, self.Src1]
@@ -122,17 +104,16 @@ class FieldsTest_Alias(unittest.TestCase):
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
         srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc([rxList0],loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc([rxList1],loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc([rxList2],loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
-        Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
+        Src3 = Survey.BaseSrc([rxList3],loc=srcLoc)
+        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3],loc=srcLoc)
         srcList = [Src0,Src1,Src2,Src3,Src4]
         survey = Survey.BaseSurvey(srcList=srcList)
-        self.D = Survey.Data(survey)
         self.F = Problem.Fields(mesh, survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',(lambda e, ind: self.F.mesh.edgeCurl * e)]})
         self.Src0 = Src0
         self.Src1 = Src1
@@ -158,7 +139,7 @@ class FieldsTest_Alias(unittest.TestCase):
 
         e = np.random.rand(F.mesh.nE,1)
         F[self.Src0, 'e'] = e
-        self.assertTrue(np.all(F[self.Src0, 'b'] == F.mesh.edgeCurl * Utils.mkvc(e)))
+        self.assertTrue(np.all(F[self.Src0, 'b'] == F.mesh.edgeCurl * e))
 
         def f():
             F[self.Src0, 'b'] = F[self.Src0, 'b']
@@ -166,7 +147,7 @@ class FieldsTest_Alias(unittest.TestCase):
 
     def test_aliasFunction(self):
         def alias(e, ind):
-            self.assertTrue(ind is self.Src0)
+            self.assertTrue(ind[0] is self.Src0)
             return self.F.mesh.edgeCurl * e
         F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
         e = np.random.rand(F.mesh.nE,1)
@@ -193,14 +174,14 @@ class FieldsTest_Time(unittest.TestCase):
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
         srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
-        Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
+        Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
+        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
         srcList = [Src0,Src1,Src2,Src3,Src4]
         survey = Survey.BaseSurvey(srcList=srcList)
         prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
@@ -249,7 +230,7 @@ class FieldsTest_Time(unittest.TestCase):
 
         b = np.random.rand(F.mesh.nF,1,nT)
         F[self.Src0, 'b', 0] = b[:,:,0]
-        self.assertTrue(np.all(F[self.Src0, 'b', 0] == b[:,0,0]))
+        self.assertTrue(np.all(F[self.Src0, 'b', 0] == Utils.mkvc(b[:,0,0],2)))
 
         phi = np.random.rand(F.mesh.nC,2,nT)
         F[[self.Src0,self.Src1], 'phi'] = phi
@@ -265,10 +246,10 @@ class FieldsTest_Time(unittest.TestCase):
         self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,nT))
         self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0,:]))
         self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1,:]))
-        self.assertTrue(np.all(F[self.Src0,'b',1] == b[:,0,1]))
-        self.assertTrue(np.all(F[self.Src1,'b',1] == b[:,1,1]))
-        self.assertTrue(np.all(F[self.Src0,'b',4] == b[:,0,4]))
-        self.assertTrue(np.all(F[self.Src1,'b',4] == b[:,1,4]))
+        self.assertTrue(np.all(F[self.Src0,'b',1] == Utils.mkvc(b[:,0,1],2)))
+        self.assertTrue(np.all(F[self.Src1,'b',1] == Utils.mkvc(b[:,1,1],2)))
+        self.assertTrue(np.all(F[self.Src0,'b',4] == Utils.mkvc(b[:,0,4],2)))
+        self.assertTrue(np.all(F[self.Src1,'b',4] == Utils.mkvc(b[:,1,4],2)))
 
 
         b = np.random.rand(F.mesh.nF, 2, nT)
@@ -295,14 +276,14 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
         srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc( [rxList0],loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc( [rxList1],loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc( [rxList2],loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
-        Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
+        Src3 = Survey.BaseSrc( [rxList3],loc=srcLoc)
+        Src4 = Survey.BaseSrc( [rxList0, rxList1, rxList2, rxList3],loc=srcLoc)
         srcList = [Src0,Src1,Src2,Src3,Src4]
         survey = Survey.BaseSurvey(srcList=srcList)
         prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
@@ -344,7 +325,7 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
         self.assertTrue(np.all(F[self.Src0, 'e', :] == e[:,0,:] ))
         self.assertTrue(np.all(F[self.Src1, 'e', :] == e[:,1,:] ))
         for t in range(nT):
-            self.assertTrue(np.all(F[self.Src1, 'e', t] == e[:,1,t] ))
+            self.assertTrue(np.all(F[self.Src1, 'e', t] == Utils.mkvc(e[:,1,t],2) ))
 
         b = np.random.rand(F.mesh.nF,nT)
         F[self.Src0, 'b',:] = b
@@ -362,7 +343,7 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
         count = [0]
         def alias(e, srcInd, timeInd):
             count[0] += 1
-            self.assertTrue(srcInd is self.Src0)
+            self.assertTrue(srcInd[0] is self.Src0)
             return self.F.mesh.edgeCurl * e
         F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
         e = np.random.rand(F.mesh.nE,1,nT)
diff --git a/SimPEG/Tests/test_LogicallyRectMesh.py b/SimPEG/Tests/test_LogicallyRectMesh.py
deleted file mode 100644
index 5d223f68..00000000
--- a/SimPEG/Tests/test_LogicallyRectMesh.py
+++ /dev/null
@@ -1,104 +0,0 @@
-import numpy as np
-import unittest
-from SimPEG.Mesh import TensorMesh, LogicallyRectMesh
-from SimPEG.Utils import ndgrid
-
-
-class BasicLRMTests(unittest.TestCase):
-
-    def setUp(self):
-        a = np.array([1, 1, 1])
-        b = np.array([1, 2])
-        c = np.array([1, 4])
-        gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h]
-        X, Y = ndgrid(gridIt([a, b]), vector=False)
-        self.TM2 = TensorMesh([a, b])
-        self.LRM2 = LogicallyRectMesh([X, Y])
-        X, Y, Z = ndgrid(gridIt([a, b, c]), vector=False)
-        self.TM3 = TensorMesh([a, b, c])
-        self.LRM3 = LogicallyRectMesh([X, Y, Z])
-
-    def test_area_3D(self):
-        test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2])
-        self.assertTrue(np.all(self.LRM3.area == test_area))
-
-    def test_vol_3D(self):
-        test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8])
-        np.testing.assert_almost_equal(self.LRM3.vol, test_vol)
-        self.assertTrue(True)  # Pass if you get past the assertion.
-
-    def test_vol_2D(self):
-        test_vol = np.array([1, 1, 1, 2, 2, 2])
-        t1 = np.all(self.LRM2.vol == test_vol)
-        self.assertTrue(t1)
-
-    def test_edge_3D(self):
-        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
-        t1 = np.all(self.LRM3.edge == test_edge)
-        self.assertTrue(t1)
-
-    def test_edge_2D(self):
-        test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2])
-        t1 = np.all(self.LRM2.edge == test_edge)
-        self.assertTrue(t1)
-
-    def test_tangents(self):
-        T = self.LRM2.tangents
-        self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM2.nEx)))
-        self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM2.nEx)))
-        self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM2.nEy)))
-        self.assertTrue(np.all(self.LRM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM2.nEy)))
-
-        T = self.LRM3.tangents
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LRM3.nEx)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LRM3.nEx)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LRM3.nEx)))
-
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LRM3.nEy)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LRM3.nEy)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LRM3.nEy)))
-
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LRM3.nEz)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LRM3.nEz)))
-        self.assertTrue(np.all(self.LRM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LRM3.nEz)))
-
-    def test_normals(self):
-        N = self.LRM2.normals
-        self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM2.nFx)))
-        self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM2.nFx)))
-        self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM2.nFy)))
-        self.assertTrue(np.all(self.LRM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM2.nFy)))
-
-        N = self.LRM3.normals
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LRM3.nFx)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LRM3.nFx)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LRM3.nFx)))
-
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LRM3.nFy)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LRM3.nFy)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LRM3.nFy)))
-
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LRM3.nFz)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LRM3.nFz)))
-        self.assertTrue(np.all(self.LRM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LRM3.nFz)))
-
-    def test_grid(self):
-        self.assertTrue(np.all(self.LRM2.gridCC == self.TM2.gridCC))
-        self.assertTrue(np.all(self.LRM2.gridN == self.TM2.gridN))
-        self.assertTrue(np.all(self.LRM2.gridFx == self.TM2.gridFx))
-        self.assertTrue(np.all(self.LRM2.gridFy == self.TM2.gridFy))
-        self.assertTrue(np.all(self.LRM2.gridEx == self.TM2.gridEx))
-        self.assertTrue(np.all(self.LRM2.gridEy == self.TM2.gridEy))
-
-        self.assertTrue(np.all(self.LRM3.gridCC == self.TM3.gridCC))
-        self.assertTrue(np.all(self.LRM3.gridN == self.TM3.gridN))
-        self.assertTrue(np.all(self.LRM3.gridFx == self.TM3.gridFx))
-        self.assertTrue(np.all(self.LRM3.gridFy == self.TM3.gridFy))
-        self.assertTrue(np.all(self.LRM3.gridFz == self.TM3.gridFz))
-        self.assertTrue(np.all(self.LRM3.gridEx == self.TM3.gridEx))
-        self.assertTrue(np.all(self.LRM3.gridEy == self.TM3.gridEy))
-        self.assertTrue(np.all(self.LRM3.gridEz == self.TM3.gridEz))
-
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/SimPEG/Tests/test_PropMaps.py b/SimPEG/Tests/test_PropMaps.py
new file mode 100644
index 00000000..b4c0eb8d
--- /dev/null
+++ b/SimPEG/Tests/test_PropMaps.py
@@ -0,0 +1,193 @@
+import unittest
+from SimPEG import *
+from scipy.constants import mu_0
+
+
+class MyPropMap(Maps.PropMap):
+    sigma = Maps.Property("Electrical Conductivity", defaultInvProp=True)
+    mu    = Maps.Property("Mu", defaultVal=mu_0)
+
+class MyReciprocalPropMap(Maps.PropMap):
+    sigma = Maps.Property("Electrical Conductivity", defaultInvProp=True, propertyLink=('rho',   Maps.ReciprocalMap))
+    rho   = Maps.Property("Electrical Resistivity",                       propertyLink=('sigma', Maps.ReciprocalMap))
+    mu    = Maps.Property("Mu", defaultVal=mu_0, propertyLink=('mui',   Maps.ReciprocalMap))
+    mui    = Maps.Property("Mu", defaultVal=1./mu_0, propertyLink=('mu',   Maps.ReciprocalMap))
+
+
+class TestPropMaps(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+    def test_setup(self):
+        expMap = Maps.ExpMap(Mesh.TensorMesh((3,)))
+        assert expMap.nP == 3
+
+        PM1 = MyPropMap(expMap)
+        PM2 = MyPropMap([('sigma', expMap)])
+        PM3 = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':slice(0,3)}})
+
+        for PM in [PM1,PM2,PM3]:
+            assert PM.defaultInvProp == 'sigma'
+            assert PM.sigmaMap is not None
+            assert PM.sigmaMap is expMap
+            assert PM.sigmaIndex == slice(0,3)
+            assert getattr(PM, 'sigma', None) is None
+            assert PM.muMap is None
+            assert PM.muIndex is None
+
+            assert 'sigma' in PM
+            assert 'mu' not in PM
+            assert 'mui' not in PM
+
+            m = PM(np.r_[1.,2,3])
+
+            assert 'sigma' in m
+            assert 'mu' not in m
+            assert 'mui' not in m
+
+            assert m.mu == mu_0
+            assert m.muModel is None
+            assert m.muMap is None
+            assert m.muDeriv is None
+
+            assert np.all(m.sigmaModel == np.r_[1.,2,3])
+            assert m.sigmaMap is expMap
+            assert np.all(m.sigma == np.exp(np.r_[1.,2,3]))
+            assert m.sigmaDeriv is not None
+
+            assert m.nP == 3
+
+    def test_slices(self):
+        expMap = Maps.ExpMap(Mesh.TensorMesh((3,)))
+        PM = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':[2,1,0]}})
+        assert PM.sigmaIndex == [2,1,0]
+        m = PM(np.r_[1.,2,3])
+        assert np.all(m.sigmaModel == np.r_[3,2,1])
+        assert np.all(m.sigma == np.exp(np.r_[3,2,1]))
+
+    def test_multiMap(self):
+        m = Mesh.TensorMesh((3,))
+        expMap = Maps.ExpMap(m)
+        iMap = Maps.IdentityMap(m)
+        PM = MyPropMap([('sigma', expMap), ('mu', iMap)])
+
+        pm = PM(np.r_[1.,2,3,4,5,6])
+
+        assert pm.nP == 6
+
+        assert 'sigma' in PM
+        assert 'mu' in PM
+        assert 'mui' not in PM
+
+        assert 'sigma' in pm
+        assert 'mu' in pm
+        assert 'mui' not in pm
+
+        assert np.all(pm.sigmaModel == [1.,2,3])
+        assert np.all(pm.sigma == np.exp([1.,2,3]))
+        assert np.all(pm.muModel == [4.,5,6])
+        assert np.all(pm.mu == [4.,5,6])
+
+
+    def test_multiMapCompressed(self):
+        m = Mesh.TensorMesh((3,))
+        expMap = Maps.ExpMap(m)
+        iMap = Maps.IdentityMap(m)
+        PM = MyPropMap({'maps':[('sigma', expMap), ('mu', iMap)],'slices':{'mu':[0,1,2]}})
+
+        pm = PM(np.r_[1,2.,3])
+
+        assert pm.nP == 3
+
+        assert 'sigma' in PM
+        assert 'mu' in PM
+        assert 'mui' not in PM
+
+        assert 'sigma' in pm
+        assert 'mu' in pm
+        assert 'mui' not in pm
+
+        assert np.all(pm.sigmaModel == [1,2,3])
+        assert np.all(pm.sigma == np.exp([1,2,3]))
+        assert np.all(pm.muModel == [1,2,3])
+        assert np.all(pm.mu == [1,2,3])
+
+    def test_Projections(self):
+        m = Mesh.TensorMesh((3,))
+        iMap = Maps.IdentityMap(m)
+        PM = MyReciprocalPropMap([('sigma', iMap)])
+        v = np.r_[1,2.,3]
+        pm = PM(v)
+
+        assert pm.sigmaProj is not None
+        assert pm.rhoProj   is None
+        assert pm.muProj    is None
+        assert pm.muiProj   is None
+
+        assert np.all(pm.sigmaProj * v == pm.sigmaModel)
+
+    def test_Links(self):
+        m = Mesh.TensorMesh((3,))
+        expMap = Maps.ExpMap(m)
+        iMap = Maps.IdentityMap(m)
+        PM = MyReciprocalPropMap([('sigma', iMap)])
+        pm = PM(np.r_[1,2.,3])
+        # print pm.sigma
+        # print pm.sigmaMap
+        assert np.all(pm.sigma == [1,2,3])
+        assert np.all(pm.rho == 1./np.r_[1,2,3])
+        assert pm.sigmaMap is iMap
+        assert pm.rhoMap is None
+        assert pm.sigmaDeriv is not None
+        assert pm.rhoDeriv is not None
+
+        assert 'sigma' in PM
+        assert 'rho' not in PM
+        assert 'mu' not in PM
+        assert 'mui' not in PM
+
+
+        assert 'sigma' in pm
+        assert 'rho' not in pm
+        assert 'mu' not in pm
+        assert 'mui' not in pm
+
+        assert pm.mu == mu_0
+        assert pm.mui == 1.0/mu_0
+        assert pm.muMap is None
+        assert pm.muDeriv is None
+        assert pm.muiMap is None
+        assert pm.muiDeriv is None
+
+        PM = MyReciprocalPropMap([('rho', iMap)])
+        pm = PM(np.r_[1,2.,3])
+        # print pm.sigma
+        # print pm.sigmaMap
+        assert np.all(pm.sigma == 1./np.r_[1,2,3])
+        assert np.all(pm.rho == [1,2,3])
+        assert pm.sigmaMap is None
+        assert pm.rhoMap is iMap
+        assert pm.sigmaDeriv is not None
+        assert pm.rhoDeriv is not None
+
+        assert 'sigma' not in PM
+        assert 'rho' in PM
+        assert 'mu' not in PM
+        assert 'mui' not in PM
+
+
+        assert 'sigma' not in pm
+        assert 'rho' in pm
+        assert 'mu' not in pm
+        assert 'mui' not in pm
+
+        self.assertRaises(AssertionError, MyReciprocalPropMap, [('rho', iMap), ('sigma', iMap)])
+        self.assertRaises(AssertionError, MyReciprocalPropMap, [('sigma', iMap), ('rho', iMap)])
+
+        MyReciprocalPropMap([('sigma', iMap), ('mu', iMap)]) # This should be fine
+
+if __name__ == '__main__':
+    unittest.main()
+
+
diff --git a/SimPEG/Tests/test_SurveyAndData.py b/SimPEG/Tests/test_SurveyAndData.py
new file mode 100644
index 00000000..f02f14e8
--- /dev/null
+++ b/SimPEG/Tests/test_SurveyAndData.py
@@ -0,0 +1,53 @@
+import unittest
+from SimPEG import *
+
+class TestData(unittest.TestCase):
+
+    def setUp(self):
+        mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
+        x = np.linspace(5,10,3)
+        XYZ = Utils.ndgrid(x,x,np.r_[0.])
+        srcLoc = np.r_[0,0,0.]
+        rxList0 = Survey.BaseRx(XYZ, 'exi')
+        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
+        rxList1 = Survey.BaseRx(XYZ, 'bxi')
+        Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
+        rxList2 = Survey.BaseRx(XYZ, 'bxi')
+        Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
+        rxList3 = Survey.BaseRx(XYZ, 'bxi')
+        Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
+        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
+        srcList = [Src0,Src1,Src2,Src3,Src4]
+        survey = Survey.BaseSurvey(srcList=srcList)
+        self.D = Survey.Data(survey)
+
+    def test_data(self):
+        V = []
+        for src in self.D.survey.srcList:
+            for rx in src.rxList:
+                v = np.random.rand(rx.nD)
+                V += [v]
+                self.D[src, rx] = v
+                self.assertTrue(np.all(v == self.D[src, rx]))
+        V = np.concatenate(V)
+        self.assertTrue(np.all(V == Utils.mkvc(self.D)))
+
+        D2 = Survey.Data(self.D.survey, V)
+        self.assertTrue(np.all(Utils.mkvc(D2) == Utils.mkvc(self.D)))
+
+    def test_uniqueSrcs(self):
+        srcs = self.D.survey.srcList
+        srcs += [srcs[0]]
+        self.assertRaises(AssertionError, Survey.BaseSurvey, srcList=srcs)
+
+    def test_sourceIndex(self):
+        survey = self.D.survey
+        srcs = survey.srcList
+        assert survey.getSourceIndex([srcs[1],srcs[0]]) == [1,0]
+        assert survey.getSourceIndex([srcs[1],srcs[2],srcs[2]]) == [1,2,2]
+        SrcNotThere = Survey.BaseSrc(srcs[0].rxList, loc=np.r_[0,0,0])
+        self.assertRaises(KeyError, survey.getSourceIndex, [SrcNotThere])
+        self.assertRaises(KeyError, survey.getSourceIndex, [srcs[1],srcs[2],SrcNotThere])
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/SimPEG/Tests/test_innerProduct.py b/SimPEG/Tests/test_innerProduct.py
index bbaf4b14..0a5ff809 100644
--- a/SimPEG/Tests/test_innerProduct.py
+++ b/SimPEG/Tests/test_innerProduct.py
@@ -7,7 +7,7 @@ from SimPEG import Utils
 class TestInnerProducts(OrderTest):
     """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
 
-    meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM']
+    meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
     meshDimension = 3
     meshSizes = [16, 32]
 
@@ -154,7 +154,7 @@ class TestInnerProducts(OrderTest):
 class TestInnerProducts2D(OrderTest):
     """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
 
-    meshTypes = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM']
+    meshTypes = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
     meshDimension = 2
     meshSizes = [4, 8, 16, 32, 64, 128]
 
diff --git a/SimPEG/Tests/test_innerProductDerivs.py b/SimPEG/Tests/test_innerProductDerivs.py
index 4b7d63cd..6eb561c1 100644
--- a/SimPEG/Tests/test_innerProductDerivs.py
+++ b/SimPEG/Tests/test_innerProductDerivs.py
@@ -7,9 +7,9 @@ from TestUtils import checkDerivative
 class TestInnerProductsDerivs(unittest.TestCase):
 
     def doTestFace(self, h, rep, fast, meshType, invProp=False, invMat=False):
-        if meshType == 'LRM':
+        if meshType == 'Curv':
             hRect = Utils.exampleLrmGrid(h,'rotate')
-            mesh = Mesh.LogicallyRectMesh(hRect)
+            mesh = Mesh.CurvilinearMesh(hRect)
         elif meshType == 'Tree':
             mesh = Mesh.TreeMesh(h)
         elif meshType == 'Tensor':
@@ -24,9 +24,9 @@ class TestInnerProductsDerivs(unittest.TestCase):
         return checkDerivative(fun, sig, num=5, plotIt=False)
 
     def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False):
-        if meshType == 'LRM':
+        if meshType == 'Curv':
             hRect = Utils.exampleLrmGrid(h,'rotate')
-            mesh = Mesh.LogicallyRectMesh(hRect)
+            mesh = Mesh.CurvilinearMesh(hRect)
         elif meshType == 'Tree':
             mesh = Mesh.TreeMesh(h)
         elif meshType == 'Tensor':
@@ -137,65 +137,65 @@ class TestInnerProductsDerivs(unittest.TestCase):
 
 
 
-    def test_FaceIP_2D_float_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],0, False, 'LRM'))
-    def test_FaceIP_3D_float_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'LRM'))
-    def test_FaceIP_2D_isotropic_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],1, False, 'LRM'))
-    def test_FaceIP_3D_isotropic_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'LRM'))
-    def test_FaceIP_2D_anisotropic_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],2, False, 'LRM'))
-    def test_FaceIP_3D_anisotropic_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'LRM'))
-    def test_FaceIP_2D_tensor_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],3, False, 'LRM'))
-    def test_FaceIP_3D_tensor_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'LRM'))
+    def test_FaceIP_2D_float_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],0, False, 'Curv'))
+    def test_FaceIP_3D_float_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Curv'))
+    def test_FaceIP_2D_isotropic_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],1, False, 'Curv'))
+    def test_FaceIP_3D_isotropic_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Curv'))
+    def test_FaceIP_2D_anisotropic_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],2, False, 'Curv'))
+    def test_FaceIP_3D_anisotropic_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Curv'))
+    def test_FaceIP_2D_tensor_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],3, False, 'Curv'))
+    def test_FaceIP_3D_tensor_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Curv'))
 
-    def test_FaceIP_2D_float_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],0, True, 'LRM'))
-    def test_FaceIP_3D_float_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'LRM'))
-    def test_FaceIP_2D_isotropic_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],1, True, 'LRM'))
-    def test_FaceIP_3D_isotropic_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'LRM'))
-    def test_FaceIP_2D_anisotropic_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4],2, True, 'LRM'))
-    def test_FaceIP_3D_anisotropic_fast_LRM(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'LRM'))
+    def test_FaceIP_2D_float_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],0, True, 'Curv'))
+    def test_FaceIP_3D_float_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Curv'))
+    def test_FaceIP_2D_isotropic_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],1, True, 'Curv'))
+    def test_FaceIP_3D_isotropic_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Curv'))
+    def test_FaceIP_2D_anisotropic_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4],2, True, 'Curv'))
+    def test_FaceIP_3D_anisotropic_fast_Curv(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Curv'))
 
-    def test_EdgeIP_2D_float_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],0, False, 'LRM'))
-    def test_EdgeIP_3D_float_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'LRM'))
-    def test_EdgeIP_2D_isotropic_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],1, False, 'LRM'))
-    def test_EdgeIP_3D_isotropic_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'LRM'))
-    def test_EdgeIP_2D_anisotropic_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],2, False, 'LRM'))
-    def test_EdgeIP_3D_anisotropic_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'LRM'))
-    def test_EdgeIP_2D_tensor_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],3, False, 'LRM'))
-    def test_EdgeIP_3D_tensor_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'LRM'))
+    def test_EdgeIP_2D_float_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],0, False, 'Curv'))
+    def test_EdgeIP_3D_float_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Curv'))
+    def test_EdgeIP_2D_isotropic_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],1, False, 'Curv'))
+    def test_EdgeIP_3D_isotropic_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Curv'))
+    def test_EdgeIP_2D_anisotropic_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],2, False, 'Curv'))
+    def test_EdgeIP_3D_anisotropic_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Curv'))
+    def test_EdgeIP_2D_tensor_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],3, False, 'Curv'))
+    def test_EdgeIP_3D_tensor_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Curv'))
 
-    def test_EdgeIP_2D_float_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],0, True, 'LRM'))
-    def test_EdgeIP_3D_float_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'LRM'))
-    def test_EdgeIP_2D_isotropic_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],1, True, 'LRM'))
-    def test_EdgeIP_3D_isotropic_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'LRM'))
-    def test_EdgeIP_2D_anisotropic_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4],2, True, 'LRM'))
-    def test_EdgeIP_3D_anisotropic_fast_LRM(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'LRM'))
+    def test_EdgeIP_2D_float_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],0, True, 'Curv'))
+    def test_EdgeIP_3D_float_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Curv'))
+    def test_EdgeIP_2D_isotropic_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],1, True, 'Curv'))
+    def test_EdgeIP_3D_isotropic_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Curv'))
+    def test_EdgeIP_2D_anisotropic_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4],2, True, 'Curv'))
+    def test_EdgeIP_3D_anisotropic_fast_Curv(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Curv'))
 
 
 
diff --git a/SimPEG/Tests/test_maps.py b/SimPEG/Tests/test_maps.py
index 9f8fef3c..a80bc40b 100644
--- a/SimPEG/Tests/test_maps.py
+++ b/SimPEG/Tests/test_maps.py
@@ -30,8 +30,8 @@ class MapTests(unittest.TestCase):
             self.assertTrue(maps.test())
 
 
-    def test_transforms_logMap(self):
-        # Note that log maps can be kinda finicky, so we are being explicit about the random seed.
+    def test_transforms_logMap_reciprocalMap(self):
+        # Note that log/reciprocal maps can be kinda finicky, so we are being explicit about the random seed.
         v2 = np.r_[ 0.40077291, 0.14410044, 0.58452314, 0.96323738, 0.01198519, 0.79754415]
         dv2 = np.r_[ 0.80653921, 0.13132446, 0.4901117, 0.03358737, 0.65473762, 0.44252488]
         v3 = np.r_[ 0.96084865, 0.34385186, 0.39430044, 0.81671285, 0.65929109, 0.2235217, 0.87897526, 0.5784033, 0.96876393, 0.63535864, 0.84130763, 0.22123854]
@@ -41,6 +41,11 @@ class MapTests(unittest.TestCase):
         maps = Maps.LogMap(self.mesh3)
         self.assertTrue(maps.test(v3, dx=dv3))
 
+        maps = Maps.ReciprocalMap(self.mesh2)
+        self.assertTrue(maps.test(v2, dx=dv2))
+        maps = Maps.ReciprocalMap(self.mesh3)
+        self.assertTrue(maps.test(v3, dx=dv3))
+
     def test_Mesh2MeshMap(self):
         maps = Maps.Mesh2Mesh([self.mesh22, self.mesh2])
         self.assertTrue(maps.test())
diff --git a/SimPEG/Tests/test_operators.py b/SimPEG/Tests/test_operators.py
index 9a465fcd..416689bc 100644
--- a/SimPEG/Tests/test_operators.py
+++ b/SimPEG/Tests/test_operators.py
@@ -4,7 +4,7 @@ from TestUtils import OrderTest
 import matplotlib.pyplot as plt
 
 #TODO: 'randomTensorMesh'
-MESHTYPES = ['uniformTensorMesh', 'uniformLRM', 'rotateLRM']
+MESHTYPES = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
 call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
 call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
 cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
@@ -38,7 +38,7 @@ class TestCurl(OrderTest):
         curlE_ana = self.M.projectFaceVector(Fc)
 
         curlE = self.M.edgeCurl.dot(E)
-        if self._meshType == 'rotateLRM':
+        if self._meshType == 'rotateCurv':
             # Really it is the integration we should be caring about:
             # So, let us look at the l2 norm.
             err = np.linalg.norm(self.M.area*(curlE - curlE_ana), 2)
@@ -229,7 +229,7 @@ class TestFaceDiv3D(OrderTest):
         divF = self.M.faceDiv.dot(F)
         divF_ana = call3(sol, self.M.gridCC)
 
-        if self._meshType == 'rotateLRM':
+        if self._meshType == 'rotateCurv':
             # Really it is the integration we should be caring about:
             # So, let us look at the l2 norm.
             err = np.linalg.norm(self.M.vol*(divF-divF_ana), 2)
diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py
index e67267cd..6b88f1d2 100644
--- a/SimPEG/Utils/__init__.py
+++ b/SimPEG/Utils/__init__.py
@@ -1,7 +1,7 @@
 from matutils import *
 from codeutils import *
 from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile
-from lrmutils import volTetra, faceInfo, indexCube
+from curvutils import volTetra, faceInfo, indexCube
 from interputils import interpmat
 from ipythonutils import easyAnimate as animate
 from CounterUtils import *
diff --git a/SimPEG/Utils/codeutils.py b/SimPEG/Utils/codeutils.py
index 1399c1b5..0ba57b2e 100644
--- a/SimPEG/Utils/codeutils.py
+++ b/SimPEG/Utils/codeutils.py
@@ -58,9 +58,11 @@ def hook(obj, method, name=None, overwrite=False, silent=False):
         print 'Method '+name+' was not overwritten.'
 
 
-def setKwargs(obj, **kwargs):
+def setKwargs(obj, ignore=[],  **kwargs):
     """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
     for attr in kwargs:
+        if attr in ignore:
+            continue
         if hasattr(obj, attr):
             setattr(obj, attr, kwargs[attr])
         else:
diff --git a/SimPEG/Utils/lrmutils.py b/SimPEG/Utils/curvutils.py
similarity index 100%
rename from SimPEG/Utils/lrmutils.py
rename to SimPEG/Utils/curvutils.py
diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py
index 475da3cc..f88f5af2 100644
--- a/SimPEG/__init__.py
+++ b/SimPEG/__init__.py
@@ -16,7 +16,7 @@ import Inversion
 import Parallel
 import Tests
 
-__version__   = '0.1.1'
+__version__   = '0.1.3'
 __author__    = 'Rowan Cockett'
 __license__   = 'MIT'
 __copyright__ = 'Copyright 2014 Rowan Cockett'
diff --git a/docs/SimPEGFrameworkRevised.png b/docs/SimPEGFrameworkRevised.png
new file mode 100644
index 00000000..18ce7a30
Binary files /dev/null and b/docs/SimPEGFrameworkRevised.png differ
diff --git a/docs/api_InnerProducts.rst b/docs/api_InnerProducts.rst
index c3fbf4d6..bc9426c6 100644
--- a/docs/api_InnerProducts.rst
+++ b/docs/api_InnerProducts.rst
@@ -75,7 +75,7 @@ We multiply by  square-root of volume on each side of the tensor conductivity to
 
 .. math::
     \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\text{TENSOR} \\
-    \mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{LRM}
+    \mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{Curv}
 
 Here the \\\(i\\\) index refers to where we choose to approximate this integral, as discussed in the note above.
 We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix \\\( \\mathbf{Q}_{(i)} \\\) to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like:
@@ -114,7 +114,7 @@ Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of
 
 .. math::
 
-    \mathbf{P}_{(i)} =  \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{LRM only}} \mathbf{Q}_{(i)}
+    \mathbf{P}_{(i)} =  \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{Curv only}} \mathbf{Q}_{(i)}
 
 .. note::
 
diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst
index c698cf67..7bf398b1 100644
--- a/docs/api_Mesh.rst
+++ b/docs/api_Mesh.rst
@@ -29,7 +29,7 @@ the implementations.
     tM = Mesh.TensorMesh(sz)
     qM = Mesh.TreeMesh(sz)
     qM.refine(lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.3 else 0)
-    rM = Mesh.LogicallyRectMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
+    rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
 
     fig, axes = plt.subplots(1,3,figsize=(14,5))
     opts = {}
@@ -38,7 +38,7 @@ the implementations.
     qM.plotGrid(ax=axes[1], **opts)
     axes[1].set_title('TreeMesh')
     rM.plotGrid(ax=axes[2], **opts)
-    axes[2].set_title('LogicallyRectMesh')
+    axes[2].set_title('CurvilinearMesh')
     plt.show()
 
 
diff --git a/docs/api_MeshCode.rst b/docs/api_MeshCode.rst
index 13ceffba..a3813c13 100644
--- a/docs/api_MeshCode.rst
+++ b/docs/api_MeshCode.rst
@@ -18,10 +18,10 @@ Tree Mesh
     :undoc-members:
 
 
-Logically Rectangular Mesh
-==========================
+Curvilinear Mesh
+================
 
-.. automodule:: SimPEG.Mesh.LogicallyRectMesh
+.. automodule:: SimPEG.Mesh.Curvilinear
     :show-inheritance:
     :members:
     :undoc-members:
diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst
index 7903cdc4..042aef67 100644
--- a/docs/api_Utils.rst
+++ b/docs/api_Utils.rst
@@ -23,10 +23,10 @@ Solver Utilities
     :members:
     :undoc-members:
 
-LRM Utilities
+Curv Utilities
 =============
 
-.. automodule:: SimPEG.Utils.lrmutils
+.. automodule:: SimPEG.Utils.curvutils
     :members:
     :undoc-members:
 
diff --git a/docs/conf.py b/docs/conf.py
index 161c663e..3bb33621 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers'
 # built documents.
 #
 # The short X.Y version.
-version = '0.1.1'
+version = '0.1.3'
 # The full version, including alpha/beta/rc tags.
-release = '0.1.1'
+release = '0.1.3'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/setup.py b/setup.py
index d1ad32e5..31c07dd5 100644
--- a/setup.py
+++ b/setup.py
@@ -33,7 +33,7 @@ with open("README.rst") as f:
 
 setup(
     name = "SimPEG",
-    version = "0.1.1",
+    version = "0.1.3",
     packages = find_packages(),
     install_requires = ['numpy>=1.7',
                         'scipy>=0.13'
@@ -44,7 +44,7 @@ setup(
     long_description = LONG_DESCRIPTION,
     license = "MIT",
     keywords = "geophysics inverse problem",
-    url = "http://simpeg.3ptscience.com/",
+    url = "http://simpeg.xyz/",
     download_url = "http://github.com/simpeg/simpeg",
     classifiers=CLASSIFIERS,
     platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"],