Index: fipy/meshes/numMesh/mesh.py
===================================================================
--- fipy/meshes/numMesh/mesh.py	(revision 2704)
+++ fipy/meshes/numMesh/mesh.py	(working copy)
@@ -48,6 +48,8 @@
 
 from fipy.tools.dimensions.physicalField import PhysicalField
 
+import numpy
+
 class MeshAdditionError(Exception):
     pass
     
@@ -202,55 +204,21 @@
         ## check dimensions
         if(self.vertexCoords.shape[0] != other.vertexCoords.shape[0]):
             raise MeshAdditionError, "Dimensions do not match"
-        ## compute vertex correlates
-        vertexCorrelates = {}
-        for i in range(selfNumVertices):
-            for j in range(otherNumVertices):
-                diff = self.vertexCoords[...,i] - other.vertexCoords[...,j]
-                diff = numerix.array(diff)
-                if (sum(diff ** 2) < smallNumber):
-                    vertexCorrelates[j] = i
+
+        vertexCorrelates = self._getVertexCorrelates(other, smallNumber)
+
         if (vertexCorrelates == {}):
             raise MeshAdditionError, "Vertices are not aligned"
 
-        
-         
-        ## compute face correlates
-        faceCorrelates = {}
-        for i in range(otherNumFaces):
-##          Seems to be overwriting other.faceVertexIDs with new numpy
-##            currFace = other.faceVertexIDs[i]
-##            currFace = other.faceVertexIDs[...,i].copy()
-##          Changed this again as numpy 1.0.4 seems to have no copy method for
-##          masked arrays.
-            try:
-                currFace = other.faceVertexIDs[...,i].copy()
-            except:
-                currFace = MA.array(other.faceVertexIDs[...,i], mask=MA.getmask(other.faceVertexIDs[...,i]))
+        faceCorrelates = self._getFaceCorrelates(other, vertexCorrelates)
 
-            keepGoing = 1
-            currIndex = 0 
-            for item in currFace:
-                if(vertexCorrelates.has_key(item)):
-                    currFace[currIndex] = vertexCorrelates[item]
-                    currIndex = currIndex + 1
-                else:
-                    keepGoing = 0
-            if(keepGoing == 1):
-                for j in range(selfNumFaces):
-                    if (self._equalExceptOrder(currFace, self.faceVertexIDs[...,j])):
-                        faceCorrelates[i] = j
         if(faceCorrelates == {}):
             raise MeshAdditionError, "Faces are not aligned"
         
-        faceIndicesToAdd = ()
-        for i in range(otherNumFaces):
-            if(not faceCorrelates.has_key(i)):
-                faceIndicesToAdd = faceIndicesToAdd + (i,)
-        vertexIndicesToAdd = ()
-        for i in range(otherNumVertices):
-            if(not vertexCorrelates.has_key(i)):
-                vertexIndicesToAdd = vertexIndicesToAdd + (i,)
+        faceIndicesToAdd = [i for i in range(otherNumFaces) 
+                              if i not in faceCorrelates]
+        vertexIndicesToAdd = [i for i in range(otherNumVertices)
+                                if i not in vertexCorrelates]
 
         ##compute the full face and vertex correlation list
         a = selfNumFaces
@@ -288,24 +256,58 @@
             'faceVertexIDs': numerix.concatenate((self.faceVertexIDs, facesToAdd), axis=1), 
             'cellFaceIDs': MA.concatenate((self.cellFaceIDs, cellsToAdd), axis=1)
             }
+            
+    def _getExteriorVertexIds(self):
+        """Return vertex ids on the exterior of the domain"""
+        fv_ids = self._getFaceVertexIDs()
+        ext = self.getExteriorFaces().getValue()
+        ext_vert_ids = numpy.unique(numpy.hstack([v[ext] for v in fv_ids]))
+        return ext_vert_ids
+            
+    def _getVertexCorrelates(self, other, smallNumber):
+        vertexCorrelates = {}
+        vertexCoords = self.getVertexCoords()
+        ext_ids = self._getExteriorVertexIds()
+        # only check exterior vertices for correlation with other
+        vertexCoords = vertexCoords[:, ext_ids]
+        for i, selfCoords in zip(ext_ids, vertexCoords.transpose()):
+            # shape selfCoords into column vector
+            selfCoords = numpy.expand_dims(selfCoords, axis=1)
+            # compute distance between selfCoords and all other.vertexCoords
+            distance = (selfCoords - other.vertexCoords)**2
+            near = numpy.where(numerix.sum(distance, axis=0) < smallNumber)
+            # `where` returns a tuple with single array. Loop over the array
+            for j in near[0]: # this should probably be: if j in near[0]
+                vertexCorrelates[j] = i
+        return vertexCorrelates
+        
+    def _getFaceCorrelates(self, other, vertexCorrelates):
+        faceCorrelates = {}
+        # sort order of vertices on *each face* (axis=0) before matching
+        sortedFaceVertexIDs = numpy.sort(self.faceVertexIDs, axis=0)
+        # MAs don't have copy method in numpy 1.0.4; use copy() if available
+        try:
+            currFace = other.faceVertexIDs[..., 0].copy()
+            copy_ma = lambda masked: masked.copy()
+        except:
+            copy_ma = lambda masked: MA.array(masked, mask=MA.getmask(masked))
+            
+        for i, currFace in enumerate(other.faceVertexIDs.transpose()):
+            # Seems to be overwriting other.faceVertexIDs with new numpy
+            currFace = copy_ma(currFace)
+            insertFace = [vertexCorrelates[item] for item in currFace 
+                                                 if item in vertexCorrelates]
+            currFace[0:len(insertFace)] = insertFace
+            if len(insertFace) == len(currFace):
+                # arrange in sorted column to match sortedFaceVertexIDs
+                column = numpy.expand_dims(numpy.sort(currFace), axis=1)
+                # bool array where all vertices surrounding a face match
+                match = numpy.all(column == sortedFaceVertexIDs, axis=0)
+                j = numpy.where(match)[0]
+                if len(j) == 1:
+                    faceCorrelates[i] = int(j)
+        return faceCorrelates
 
-    def _equalExceptOrder(self, first, second):
-        """Determines if two lists contain the same set of elements, although they may be in different orders. Does not work if one list contains duplicates of an element.
-        """
-        res = 0
- 
-        if (len(first) == len(second)):
-            res = 1
-        for i in first:
-            isthisin = 0
-            for j in second:
-                if (i == j):
-                    isthisin = 1
-            if(isthisin == 0):
-                res = 0
-        return res
-    
-
     def _translate(self, vector):
         newCoords = self.vertexCoords + vector
         newmesh = Mesh(newCoords, numerix.array(self.faceVertexIDs), numerix.array(self.cellFaceIDs))