diff --git a/analysis/analysis.cpp b/analysis/analysis.cpp
index bc6f40f2..f28827e9 100644
--- a/analysis/analysis.cpp
+++ b/analysis/analysis.cpp
@@ -460,6 +460,7 @@ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& r
 /******************************************************************
 * Compute the mapping of blob ids between timesteps               *
 ******************************************************************/
+typedef std::map<BlobIDType,std::map<BlobIDType,int64_t> > map_type;
 template<class TYPE> inline MPI_Datatype getMPIType();
 template<> inline MPI_Datatype getMPIType<int32_t>() { return MPI_INT; }
 template<> inline MPI_Datatype getMPIType<int64_t>() { 
@@ -485,66 +486,87 @@ void gatherSet( std::set<TYPE>& set, MPI_Comm comm )
     for (size_t i=0; i<recv_data.size(); i++)
         set.insert(recv_data[i]);
 }
-template<class TYPE>
-void gatherSrcIDMap( std::map<TYPE,std::set<TYPE> >& src_map, MPI_Comm comm )
+void gatherSrcIDMap( map_type& src_map, MPI_Comm comm )
 {
     int nprocs = comm_size(comm);
-    MPI_Datatype type = getMPIType<TYPE>();
-    std::vector<TYPE> send_data;
-    typename std::map<TYPE,std::set<TYPE> >::const_iterator it;
-    for (it=src_map.begin(); it!=src_map.end(); ++it) {
+    MPI_Datatype type = getMPIType<int64_t>();
+    std::vector<int64_t> send_data;
+    for (map_type::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) {
         int id = it->first;
-        const std::set<TYPE>& src_ids = it->second;
+        const std::map<BlobIDType,int64_t>& src_ids = it->second;
         send_data.push_back(id);
         send_data.push_back(src_ids.size());
-        typename std::set<TYPE>::const_iterator it2;
-        for (it2=src_ids.begin(); it2!=src_ids.end(); ++it2)
-            send_data.push_back(*it2);
+        typename std::map<BlobIDType,int64_t>::const_iterator it2;
+        for (it2=src_ids.begin(); it2!=src_ids.end(); ++it2) {
+            send_data.push_back(it2->first);
+            send_data.push_back(it2->second);
+        }
     }
     int send_count = send_data.size();
     std::vector<int> recv_count(nprocs,0), recv_disp(nprocs,0);
     MPI_Allgather(&send_count,1,MPI_INT,getPtr(recv_count),1,MPI_INT,comm);
     for (int i=1; i<nprocs; i++)
         recv_disp[i] = recv_disp[i-1] + recv_count[i-1];
-    std::vector<TYPE> recv_data(recv_disp[nprocs-1]+recv_count[nprocs-1]);
+    std::vector<int64_t> recv_data(recv_disp[nprocs-1]+recv_count[nprocs-1]);
     MPI_Allgatherv(getPtr(send_data),send_count,type,
         getPtr(recv_data),getPtr(recv_count),getPtr(recv_disp),type,comm);
     size_t i=0;
+    src_map.clear();
     while ( i < recv_data.size() ) {
-        int id = recv_data[i];
-        int count = recv_data[i+1];
+        BlobIDType id = recv_data[i];
+        size_t count = recv_data[i+1];
         i += 2;
-        std::set<TYPE>& src_ids = src_map[id];
-        for (int j=0; j<count; j++,i++)
-            src_ids.insert(recv_data[i]);
+        std::map<BlobIDType,int64_t>& src_ids = src_map[id];
+        for (size_t j=0; j<count; j++,i+=2) {
+            std::map<BlobIDType,int64_t>::iterator it = src_ids.find(recv_data[i]);
+            if ( it == src_ids.end() )
+                src_ids.insert(std::pair<BlobIDType,int64_t>(recv_data[i],recv_data[i+1]));
+            else
+                it->second += recv_data[i+1];
+        }
     }
 }
-void addSrcDstIDs( BlobIDType src_id, std::map<BlobIDType,std::set<BlobIDType> >& src_map, 
-    std::map<BlobIDType,std::set<BlobIDType> >& dst_map, std::set<BlobIDType>& src, std::set<BlobIDType>& dst )
+void addSrcDstIDs( BlobIDType src_id, map_type& src_map, map_type& dst_map, 
+    std::set<BlobIDType>& src, std::set<BlobIDType>& dst )
 {
     src.insert(src_id);
-    const std::set<BlobIDType>& dst_ids = dst_map[src_id];
-    for (std::set<BlobIDType>::const_iterator it=dst_ids.begin(); it!=dst_ids.end(); ++it) {
-        if ( dst.find(*it)==dst.end() )
-            addSrcDstIDs(*it,dst_map,src_map,dst,src);
+    const std::map<BlobIDType,int64_t>& dst_ids = dst_map[src_id];
+    for (std::map<BlobIDType,int64_t>::const_iterator it=dst_ids.begin(); it!=dst_ids.end(); ++it) {
+        if ( dst.find(it->first)==dst.end() )
+            addSrcDstIDs(it->first,dst_map,src_map,dst,src);
     }
 }
-ID_map_struct computeIDMap( const BlobIDArray& ID1, const BlobIDArray& ID2, MPI_Comm comm )
+ID_map_struct computeIDMap( int nx, int ny, int nz, 
+    const BlobIDArray& ID1, const BlobIDArray& ID2, MPI_Comm comm )
 {
     ASSERT(ID1.size()==ID2.size());
     PROFILE_START("computeIDMap");
+    const int ngx = (ID1.size(0)-nx)/2;
+    const int ngy = (ID1.size(1)-ny)/2;
+    const int ngz = (ID1.size(2)-nz)/2;
 
     // Get a global list of all src/dst ids and the src map for each local blob
     std::set<BlobIDType> src_set, dst_set;
-    std::map<BlobIDType,std::set<BlobIDType> > src_map;   // Map of the src ids for each dst id
-    for (size_t i=0; i<ID1.length(); i++) {
-        if ( ID1(i)>=0 )
-            src_set.insert(ID1(i));
-        if ( ID2(i)>=0 )
-            dst_set.insert(ID2(i));
-        if ( ID2(i)>=0 && ID1(i)>=0 ) {
-            std::set<BlobIDType>& src_ids = src_map[ID2(i)];
-            src_ids.insert(ID1(i));
+    map_type src_map;   // Map of the src ids for each dst id
+    for (int k=ngz; k<ngz+nz; k++) {
+        for (int j=ngy; j<ngy+ny; j++) {
+            for (int i=ngx; i<ngx+nx; i++) {
+                int id1 = ID1(i,j,k);
+                int id2 = ID2(i,j,k);
+                if ( id1>=0 )
+                    src_set.insert(id1);
+                if ( id2>=0 )
+                    dst_set.insert(id2);
+                if ( id1>=0 && id2>=0 ) {
+                    std::map<BlobIDType,int64_t>& src_ids = src_map[id2];
+                    std::map<BlobIDType,int64_t>::iterator it = src_ids.find(id1);
+                    if ( it == src_ids.end() ) {
+                        src_ids.insert(std::pair<BlobIDType,int64_t>(id1,0));
+                        it = src_ids.find(id1);
+                    }
+                    it->second++;
+                }
+            }
         }
     }
     // Communicate the src/dst ids and src id map to all processors and reduce
@@ -552,13 +574,13 @@ ID_map_struct computeIDMap( const BlobIDArray& ID1, const BlobIDArray& ID2, MPI_
     gatherSet( dst_set, comm );
     gatherSrcIDMap( src_map, comm );
     // Compute the dst id map
-    std::map<BlobIDType,std::set<BlobIDType> > dst_map;   // Map of the dst ids for each src id
-    for (std::map<BlobIDType,std::set<BlobIDType> >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) {
+    map_type dst_map;   // Map of the dst ids for each src id
+    for (map_type::const_iterator it=src_map.begin(); it!=src_map.end(); ++it) {
         BlobIDType id = it->first;
-        const std::set<BlobIDType>& src_ids = it->second;
-        for (std::set<BlobIDType>::const_iterator it2=src_ids.begin(); it2!=src_ids.end(); ++it2) {
-            std::set<BlobIDType>& dst_ids = dst_map[*it2];
-            dst_ids.insert(id);
+        const std::map<BlobIDType,int64_t>& src_ids = it->second;
+        for (std::map<BlobIDType,int64_t>::const_iterator it2=src_ids.begin(); it2!=src_ids.end(); ++it2) {
+            std::map<BlobIDType,int64_t>& dst_ids = dst_map[it2->first];
+            dst_ids.insert(std::pair<BlobIDType,int64_t>(id,it2->second));
         }
     }
 
@@ -577,16 +599,16 @@ ID_map_struct computeIDMap( const BlobIDArray& ID1, const BlobIDArray& ID2, MPI_
     // Find blobs with a 1-to-1 mapping
     std::vector<BlobIDType> dst_list;
     dst_list.reserve(src_map.size());
-    for (std::map<BlobIDType,std::set<BlobIDType> >::const_iterator it=src_map.begin(); it!=src_map.end(); ++it)
+    for (map_type::const_iterator it=src_map.begin(); it!=src_map.end(); ++it)
         dst_list.push_back(it->first);
     for (size_t i=0; i<dst_list.size(); i++) {
         int dst_id = dst_list[i];
-        const std::set<BlobIDType>& src_ids = src_map[dst_id];
+        const std::map<BlobIDType,int64_t>& src_ids = src_map[dst_id];
         if ( src_ids.size()==1 ) {
-            int src_id = *src_ids.begin();
-            const std::set<BlobIDType>& dst_ids = dst_map[src_id];
+            int src_id = src_ids.begin()->first;
+            const std::map<BlobIDType,int64_t>& dst_ids = dst_map[src_id];
             if ( dst_ids.size()==1 ) {
-                ASSERT(*dst_ids.begin()==dst_id);
+                ASSERT(dst_ids.begin()->first==dst_id);
                 src_map.erase(dst_id);
                 dst_map.erase(src_id);
                 id_map.src_dst.push_back(std::pair<BlobIDType,BlobIDType>(src_id,dst_id));
@@ -598,10 +620,6 @@ ID_map_struct computeIDMap( const BlobIDArray& ID1, const BlobIDArray& ID2, MPI_
         // Get a lit of the src-dst ids
         std::set<BlobIDType> src, dst;
         addSrcDstIDs( dst_map.begin()->first, src_map, dst_map, src, dst );
-        for (std::set<BlobIDType>::const_iterator it=src.begin(); it!=src.end(); ++it)
-            dst_map.erase(*it);
-        for (std::set<BlobIDType>::const_iterator it=dst.begin(); it!=dst.end(); ++it)
-            src_map.erase(*it);
         if ( src.size()==1 ) {
             // Bubble split
             id_map.split.push_back( BlobIDSplitStruct(*src.begin(),std::vector<BlobIDType>(dst.begin(),dst.end())) );
@@ -613,6 +631,22 @@ ID_map_struct computeIDMap( const BlobIDArray& ID1, const BlobIDArray& ID2, MPI_
             id_map.merge_split.push_back( BlobIDMergeSplitStruct(
                 std::vector<BlobIDType>(src.begin(),src.end()), std::vector<BlobIDType>(dst.begin(),dst.end() ) ) );
         }
+        // Add the overlaps
+        for (std::set<BlobIDType>::const_iterator it1=src.begin(); it1!=src.end(); ++it1) {
+            const std::map<BlobIDType,int64_t>& dst_ids = dst_map[*it1];
+            for (std::set<BlobIDType>::const_iterator it2=dst.begin(); it2!=dst.end(); ++it2) {
+                std::pair<BlobIDType,BlobIDType> id(*it1,*it2);
+                int64_t overlap = 0;
+                const std::map<BlobIDType,int64_t>::const_iterator it = dst_ids.find(*it2);
+                if ( it != dst_ids.end() ) { overlap = it->second; }
+                id_map.overlap.insert(std::pair<OverlapID,int64_t>(id,overlap));
+            }
+        }
+        // Clear the mapped entries
+        for (std::set<BlobIDType>::const_iterator it=src.begin(); it!=src.end(); ++it)
+            dst_map.erase(*it);
+        for (std::set<BlobIDType>::const_iterator it=dst.begin(); it!=dst.end(); ++it)
+            src_map.erase(*it);
     }
     ASSERT(src_map.empty());
 
@@ -624,63 +658,109 @@ ID_map_struct computeIDMap( const BlobIDArray& ID1, const BlobIDArray& ID2, MPI_
 /******************************************************************
 * Renumber the ids                                                *
 ******************************************************************/
+typedef std::vector<BlobIDType> IDvec;
+inline void renumber( const std::vector<BlobIDType>& id1, const std::vector<BlobIDType>& id2,
+    const std::map<OverlapID,int64_t>& overlap, std::vector<BlobIDType>& new_ids, BlobIDType& id_max )
+{
+    if ( id2.empty() ) {
+        // No dst ids to set
+    } else if ( id1.empty() ) {
+        // No src ids
+        for (size_t i=0; i<id2.size(); i++) {
+            id_max++;
+            if ( (BlobIDType) new_ids.size() < id2[i]+1 )
+                new_ids.resize(id2[i],-1);
+            new_ids[id2[i]] = id_max;
+        }
+    } else if ( id1.size()==1 && id2.size()==1 ) {
+        // Direct src-dst mapping
+        if ( (BlobIDType) new_ids.size() < id2[0]+1 )
+            new_ids.resize(id2[0]+1,-1);
+        new_ids[id2[0]] = id1[0];
+    } else {
+        // General N to M mapping
+        // Get the overlap weights
+        Array<int64_t> cost(id1.size(),id2.size());
+        for (size_t j=0; j<id2.size(); j++) {
+            for (size_t i=0; i<id1.size(); i++) {
+                cost(i,j) = overlap.find(std::pair<BlobIDType,BlobIDType>(id1[i],id2[j]))->second;
+            }
+        }
+        // While we have not mapped all dst ids
+        while ( 1 ) {
+            size_t index = 1;
+            int64_t cost2 = -1;
+            for (size_t i=0; i<cost.length(); i++) {
+                if ( cost(i) > cost2 ) {
+                    cost2 = cost(i);
+                    index = i;
+                }
+            }
+            if ( cost2 <= 0 )
+                break;
+            // Use id1[i] for id2[j]
+            int i = index%id1.size();
+            int j = index/id1.size();
+            if ( (BlobIDType) new_ids.size() < id2[j]+1 )
+                new_ids.resize(id2[j],-1);
+            new_ids[id2[j]] = id1[i];
+            for (size_t k=0; k<id2.size(); k++)
+                cost(i,k) = -1;
+            for (size_t k=0; k<id1.size(); k++)
+                cost(k,j) = -1;
+        }
+        // No remaining src overlap with dst, create new ids for all remaining dst
+        for (size_t i=0; i<id2.size(); i++) {
+            if ( (BlobIDType) new_ids.size() < id2[i]+1 )
+                new_ids.resize(id2[i]+1,-1);
+            if ( new_ids[id2[i]] == -1 ) {
+                id_max++;
+                new_ids[id2[i]] = id_max;
+            }
+        }
+    }
+}
+inline void renumberIDs( const std::vector<BlobIDType>& new_ids, BlobIDType& id )
+{
+    id = new_ids[id];
+}
+inline void renumberIDs( const std::vector<BlobIDType>& new_ids, std::vector<BlobIDType>& ids )
+{
+    for (size_t i=0; i<ids.size(); i++)
+        ids[i] = new_ids[ids[i]];
+}
 void getNewIDs( ID_map_struct& map, BlobIDType& id_max, std::vector<BlobIDType>& new_ids )
 {
     new_ids.resize(0);
-    // Renumber the ids that map directly
-    for (size_t i=0; i<map.src_dst.size(); i++) {
-        int id1 = map.src_dst[i].second;
-        int id2 = map.src_dst[i].first;
-        map.src_dst[i].second = id2;
-        if ( new_ids.size() < static_cast<size_t>(id1+1) )
-            new_ids.resize(id1+1,-1);
-        new_ids[id1] = id2;
+    // Get the new id numbers for each map type
+    for (size_t i=0; i<map.src_dst.size(); i++)
+        renumber(IDvec(1,map.src_dst[i].first),IDvec(1,map.src_dst[i].second),map.overlap,new_ids,id_max);
+    for (size_t i=0; i<map.created.size(); i++)
+        renumber(std::vector<BlobIDType>(),IDvec(1,map.created[i]),map.overlap,new_ids,id_max);
+    for (size_t i=0; i<map.destroyed.size(); i++)
+        renumber(IDvec(1,map.destroyed[i]),std::vector<BlobIDType>(),map.overlap,new_ids,id_max);
+    for (size_t i=0; i<map.split.size(); i++)
+        renumber(IDvec(1,map.split[i].first),map.split[i].second,map.overlap,new_ids,id_max);
+    for (size_t i=0; i<map.merge.size(); i++)
+        renumber(map.merge[i].first,IDvec(1,map.merge[i].second),map.overlap,new_ids,id_max);
+    for (size_t i=0; i<map.merge_split.size(); i++)
+        renumber(map.merge_split[i].first,map.merge_split[i].second,map.overlap,new_ids,id_max);
+    // Renumber the ids in the map
+    for (size_t i=0; i<map.src_dst.size(); i++)
+        renumberIDs( new_ids, map.src_dst[i].second );
+    renumberIDs( new_ids, map.created );
+    for (size_t i=0; i<map.split.size(); i++)
+        renumberIDs( new_ids, map.split[i].second );
+    for (size_t i=0; i<map.merge.size(); i++)
+        renumberIDs( new_ids, map.merge[i].second );
+    for (size_t i=0; i<map.merge_split.size(); i++)
+        renumberIDs( new_ids, map.merge_split[i].second );
+    std::map<OverlapID,int64_t> overlap2;
+    for (std::map<OverlapID,int64_t>::const_iterator it=map.overlap.begin(); it!=map.overlap.begin(); ++it) {
+        OverlapID id = it->first;
+        renumberIDs( new_ids, id.second );
+        overlap2.insert( std::pair<OverlapID,int64_t>(id,it->second) );
     }
-    // Renumber the created blobs to create new ids
-    for (size_t i=0; i<map.created.size(); i++) {
-        int id1 = map.created[i];
-        id_max++;
-        int id2 = id_max;
-        map.created[i] = id2;
-        if ( new_ids.size() < static_cast<size_t>(id1+1) )
-            new_ids.resize(id1+1,-1);
-        new_ids[id1] = id2;
-    }
-    // Renumber the blob splits to create new ids
-    for (size_t i=0; i<map.split.size(); i++) {
-        for (size_t j=0; j<map.split[i].second.size(); j++) {
-            int id1 = map.split[i].second[j];
-            id_max++;
-            int id2 = id_max;
-            map.split[i].second[j] = id2;
-            if ( new_ids.size() < static_cast<size_t>(id1+1) )
-                new_ids.resize(id1+1,-1);
-            new_ids[id1] = id2;
-        }
-    }
-    // Renumber the blob merges to create a new id
-    for (size_t i=0; i<map.merge.size(); i++) {
-        int id1 = map.merge[i].second;
-        id_max++;
-        int id2 = id_max;
-        map.merge[i].second = id2;
-        if ( new_ids.size() < static_cast<size_t>(id1+1) )
-            new_ids.resize(id1+1,-1);
-        new_ids[id1] = id2;
-    }
-    // Renumber the blob merge/splits to create new ids
-    for (size_t i=0; i<map.merge_split.size(); i++) {
-        for (size_t j=0; j<map.merge_split[i].second.size(); j++) {
-            int id1 = map.merge_split[i].second[j];
-            id_max++;
-            int id2 = id_max;
-            map.merge_split[i].second[j] = id2;
-            if ( new_ids.size() < static_cast<size_t>(id1+1) )
-                new_ids.resize(id1+1,-1);
-            new_ids[id1] = id2;
-        }
-    }
-
 }
 void renumberIDs( const std::vector<BlobIDType>& new_ids, BlobIDArray& IDs )
 {
diff --git a/analysis/analysis.h b/analysis/analysis.h
index 2bc34128..dc73314a 100644
--- a/analysis/analysis.h
+++ b/analysis/analysis.h
@@ -93,13 +93,15 @@ void ReorderBlobIDs( BlobIDArray& ID, MPI_Comm comm );
 typedef std::pair<BlobIDType,std::vector<BlobIDType> > BlobIDSplitStruct;
 typedef std::pair<std::vector<BlobIDType>,BlobIDType> BlobIDMergeStruct;
 typedef std::pair<std::vector<BlobIDType>,std::vector<BlobIDType> > BlobIDMergeSplitStruct;
+typedef std::pair<BlobIDType,BlobIDType> OverlapID;
 struct ID_map_struct {
-    std::vector<BlobIDType> created;                   // list of new blobs that were created
-    std::vector<BlobIDType> destroyed;                 // list of blobs that disappeared
+    std::vector<BlobIDType> created;            // list of new blobs that were created
+    std::vector<BlobIDType> destroyed;          // list of blobs that disappeared
     std::vector<std::pair<BlobIDType,BlobIDType> > src_dst;   // one-one mapping of blobs (first,second timestep id)
     std::vector<BlobIDSplitStruct> split;       // list of blobs that split
     std::vector<BlobIDMergeStruct> merge;       // list of blobs that merged
     std::vector<BlobIDMergeSplitStruct> merge_split; // list of blobs that both merged and split
+    std::map<OverlapID,int64_t> overlap;         // for ids that are not a 1-1 mapping, this is a list of the overlaps <src,dst>
     //! Empty constructor
     ID_map_struct() {}
     //! Create initial map from N blobs (ordered 1:N-1)
@@ -118,7 +120,7 @@ struct ID_map_struct {
  * @param[in] ID1           The blob ids at the first timestep
  * @param[in] ID2           The blob ids at the second timestep
  */
-ID_map_struct computeIDMap( const BlobIDArray& ID1, const BlobIDArray& ID2, MPI_Comm comm );
+ID_map_struct computeIDMap( int nx, int ny, int nz, const BlobIDArray& ID1, const BlobIDArray& ID2, MPI_Comm comm );
 
 
 /*!
diff --git a/cmake/macros.cmake b/cmake/macros.cmake
index 39f1ee5d..cb068d05 100644
--- a/cmake/macros.cmake
+++ b/cmake/macros.cmake
@@ -26,6 +26,15 @@ IF ( NOT TARGET copy-${PROJ}-include )
 ENDIF()
 
 
+# Dummy use to prevent unused cmake variable warning
+MACRO( NULL_USE VAR )
+    IF ( "${${VAR}}" STREQUAL "dummy" )
+        MESSAGE( FATAL_ERROR "NULL_USE fail" )
+    ENDIF()
+ENDMACRO()
+NULL_USE( CMAKE_C_FLAGS )
+
+
 # Macro to set a global variable
 MACRO(GLOBAL_SET VARNAME)
   SET(${VARNAME} ${ARGN} CACHE INTERNAL "")
@@ -222,6 +231,7 @@ MACRO( INSTALL_${PROJ}_TARGET PACKAGE )
             CUDA_COMPILE( CUBINS ${CUDASOURCES} )
         ENDIF()
         ADD_LIBRARY( ${PACKAGE} ${LIB_TYPE} ${SOURCES} ${CUBINS} )
+        TARGET_LINK_LIBRARIES( ${PACKAGE} ${SYSTEM_LDFLAGS} )
         IF ( TARGET write_repo_version )
             ADD_DEPENDENCIES( ${PACKAGE} write_repo_version )
         ENDIF()
@@ -281,32 +291,32 @@ MACRO( IDENTIFY_COMPILER )
     IF ( CMAKE_C_COMPILER_WORKS OR CMAKE_C_COMPILER_WORKS )
         IF( CMAKE_COMPILE_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX )
             SET( USING_GCC TRUE )
-            ADD_DEFINITIONS( "-D USING_GCC" )
+            ADD_DEFINITIONS( -DUSING_GCC )
             MESSAGE("Using gcc")
         ELSEIF( MSVC OR MSVC_IDE OR MSVC60 OR MSVC70 OR MSVC71 OR MSVC80 OR CMAKE_COMPILER_2005 OR MSVC90 OR MSVC10 )
             IF( NOT ${CMAKE_SYSTEM_NAME} STREQUAL "Windows" )
                 MESSAGE( FATAL_ERROR "Using microsoft compilers on non-windows system?" )
             ENDIF()
             SET( USING_MSVC TRUE )
-            ADD_DEFINITIONS( "-D USING_MSVC" )
+            ADD_DEFINITIONS( -DUSING_MSVC )
             MESSAGE("Using Microsoft")
         ELSEIF( (${CMAKE_C_COMPILER_ID} MATCHES "Intel") OR (${CMAKE_CXX_COMPILER_ID} MATCHES "Intel") ) 
             SET(USING_ICC TRUE)
-            ADD_DEFINITIONS( "-D USING_ICC" )
+            ADD_DEFINITIONS( -DUSING_ICC )
             MESSAGE("Using icc")
         ELSEIF( (${CMAKE_C_COMPILER_ID} MATCHES "PGI") OR (${CMAKE_CXX_COMPILER_ID} MATCHES "PGI") )
             SET(USING_PGCC TRUE)
-            ADD_DEFINITIONS( "-D USING_PGCC" )
+            ADD_DEFINITIONS( -DUSING_PGCC )
             MESSAGE("Using pgCC")
         ELSEIF( (${CMAKE_C_COMPILER_ID} MATCHES "CRAY") OR (${CMAKE_CXX_COMPILER_ID} MATCHES "CRAY") OR
                 (${CMAKE_C_COMPILER_ID} MATCHES "Cray") OR (${CMAKE_CXX_COMPILER_ID} MATCHES "Cray") )
             SET(USING_CRAY TRUE)
-            ADD_DEFINITIONS( "-D USING_CRAY" )
+            ADD_DEFINITIONS( -DUSING_CRAY )
             MESSAGE("Using Cray")
         ELSEIF( (${CMAKE_C_COMPILER_ID} MATCHES "CLANG") OR (${CMAKE_CXX_COMPILER_ID} MATCHES "CLANG") OR
                 (${CMAKE_C_COMPILER_ID} MATCHES "Clang") OR (${CMAKE_CXX_COMPILER_ID} MATCHES "Clang") )
             SET(USING_CLANG TRUE)
-            ADD_DEFINITIONS( "-D USING_CLANG" )
+            ADD_DEFINITIONS( -DUSING_CLANG )
             MESSAGE("Using Clang")
         ELSE()
             SET(USING_DEFAULT TRUE)
@@ -338,6 +348,8 @@ ENDMACRO()
 MACRO( SET_WARNINGS )
   IF ( USING_GCC )
     # Add gcc specific compiler options
+    # Note: adding -Wlogical-op causes a wierd linking error on Titan using the nvcc wrapper:
+    #    /usr/bin/ld: cannot find gical-op: No such file or directory
     SET(CMAKE_C_FLAGS     "${CMAKE_C_FLAGS} -Wall -Wextra -Wcast-align -Wlogical-op -Wno-char-subscripts -Wno-unused-parameter -Wno-unused-variable") 
     SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wcast-align -Wlogical-op -Wno-char-subscripts -Wno-unused-parameter -Wno-unused-variable")
   ELSEIF ( USING_MSVC )
@@ -372,9 +384,9 @@ ENDMACRO()
 
 # Macro to add user compile flags
 MACRO( ADD_USER_FLAGS )
-    SET(CMAKE_C_FLAGS   " ${CMAKE_C_FLAGS} ${CFLAGS} ${CFLAGS_EXTRA}" )
-    SET(CMAKE_CXX_FLAGS " ${CMAKE_CXX_FLAGS} ${CXXFLAGS} ${CXXFLAGS_EXTRA}" )
-    SET(CMAKE_Fortran_FLAGS " ${CMAKE_Fortran_FLAGS} ${FFLAGS} ${FFLAGS_EXTRA}" )
+    STRING( STRIP "${CMAKE_C_FLAGS} ${CFLAGS} ${CFLAGS_EXTRA}" CMAKE_C_FLAGS )
+    STRING( STRIP "${CMAKE_CXX_FLAGS} ${CXXFLAGS} ${CXXFLAGS_EXTRA}" CMAKE_CXX_FLAGS )
+    STRING( STRIP "${CMAKE_Fortran_FLAGS} ${FFLAGS} ${FFLAGS_EXTRA}" CMAKE_Fortran_FLAGS )
 ENDMACRO()
 
 
@@ -390,44 +402,52 @@ MACRO( ADD_CXX_STD )
         MESSAGE( FATAL_ERROR "Unknown c++ standard (98,11,14,NONE)" )
     ENDIF()
     # Add the flags
-    ADD_DEFINITIONS( -DCXX_STD=${CXX_STD} )
     SET( CMAKE_CXX_STANDARD ${CXX_STD} )
     MESSAGE( "C++ standard: ${CXX_STD}" )
+    SET( CXX_STD_FLAG )
     IF ( USING_GCC )
         # GNU: -std=
         IF ( ${CXX_STD} STREQUAL "98" )
-            SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++98")
+            SET( CXX_STD_FLAG -std=c++98 )
         ELSEIF ( ${CXX_STD} STREQUAL "11" )
-            SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+            SET( CXX_STD_FLAG -std=c++11 )
         ELSEIF ( ${CXX_STD} STREQUAL "14" )
-            SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++1y")
+            SET( CXX_STD_FLAG -std=c++1y )
+        ELSE()
+            MESSAGE(FATAL_ERROR "Unknown standard")
         ENDIF()
     ELSEIF ( USING_MSVC )
         # Microsoft: Does not support this level of control
     ELSEIF ( USING_ICC )
         # ICC: -std=
-        SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++${CXX_STD}")
+        SET( CXX_STD_FLAG -std=c++${CXX_STD} )
     ELSEIF ( USING_CRAY )
         # Cray: Does not seem to support controlling the std?
     ELSEIF ( USING_PGCC )
         # PGI: -std=
         IF ( ${CXX_STD} STREQUAL "98" )
-            SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --c++0x")
+            SET( CXX_STD_FLAG --c++0x )
         ELSEIF ( ${CXX_STD} STREQUAL "11" )
-            SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --c++11")
+            SET( CXX_STD_FLAG --c++11 )
         ELSEIF ( ${CXX_STD} STREQUAL "14" )
             MESSAGE( FATAL_ERROR "C++14 features are not availible yet for PGI" )
+        ELSE()
+            MESSAGE(FATAL_ERROR "Unknown standard")
         ENDIF()
     ELSEIF ( USING_CLANG )
         # Clang: -std=
-        IF ( ( ${CXX_STD} STREQUAL "98") OR ( ${CXX_STD} STREQUAL "98" ) )
-            SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++${CXX_STD}")
+        IF ( ( ${CXX_STD} STREQUAL "98") OR ( ${CXX_STD} STREQUAL "11" ) )
+            SET( CXX_STD_FLAG -std=c++${CXX_STD} )
         ELSEIF ( ${CXX_STD} STREQUAL "14" )
-            SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++1y")
+            SET( CXX_STD_FLAG -std=c++1y )
+        ELSE()
+            MESSAGE(FATAL_ERROR "Unknown standard")
         ENDIF()
     ELSEIF ( USING_DEFAULT )
         # Default: do nothing
     ENDIF()
+    ADD_DEFINITIONS( -DCXX_STD=${CXX_STD} )
+    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CXX_STD_FLAG}")
 ENDMACRO()
 
 
@@ -437,35 +457,11 @@ MACRO( SET_COMPILER_FLAGS )
     IDENTIFY_COMPILER()
     # Set the default flags for each build type
     IF ( USING_MSVC )
-        SET(CMAKE_C_FLAGS_DEBUG       "-D_DEBUG /DEBUG /Od /EHsc /MDd /Zi /Z7" )
-        SET(CMAKE_C_FLAGS_RELEASE     "/O2 /EHsc /MD"                      )
-        SET(CMAKE_CXX_FLAGS_DEBUG     "-D_DEBUG /DEBUG /Od /EHsc /MDd /Zi /Z7" )
-        SET(CMAKE_CXX_FLAGS_RELEASE   "/O2 /EHsc /MD"                      )
-        SET(CMAKE_Fortran_FLAGS_DEBUG ""                                   )
-        SET(CMAKE_Fortran_FLAGS_RELEASE ""                                 )
+        SET(CMAKE_C_FLAGS_DEBUG   "-D_DEBUG /DEBUG /Od /EHsc /MDd /Zi /Z7" )
+        SET(CMAKE_CXX_FLAGS_DEBUG "-D_DEBUG /DEBUG /Od /EHsc /MDd /Zi /Z7" )
+        SET(CMAKE_C_FLAGS_RELEASE   "/O2 /EHsc /MD" )
+        SET(CMAKE_CXX_FLAGS_RELEASE "/O2 /EHsc /MD" )
     ELSE()
-        SET(CMAKE_C_FLAGS_DEBUG       "-g -D_DEBUG -O0" )
-        SET(CMAKE_C_FLAGS_RELEASE     "-O2"             )
-        SET(CMAKE_CXX_FLAGS_DEBUG     "-g -D_DEBUG -O0" )
-        SET(CMAKE_CXX_FLAGS_RELEASE   "-O2"             )
-        SET(CMAKE_Fortran_FLAGS_DEBUG "-g -O0"          )
-        SET(CMAKE_Fortran_FLAGS_RELEASE "-O2"           )
-    ENDIF()
-    # Set the compiler flags to use
-    IF ( ${CMAKE_BUILD_TYPE} STREQUAL "Debug" OR ${CMAKE_BUILD_TYPE} STREQUAL "DEBUG")
-        SET(CMAKE_C_FLAGS       ${CMAKE_C_FLAGS_DEBUG}       )
-        SET(CMAKE_CXX_FLAGS     ${CMAKE_CXX_FLAGS_DEBUG}     )
-        SET(CMAKE_Fortran_FLAGS ${CMAKE_Fortran_FLAGS_DEBUG} )
-    ELSEIF ( ${CMAKE_BUILD_TYPE} STREQUAL "Release" OR ${CMAKE_BUILD_TYPE} STREQUAL "RELEASE")
-        SET(CMAKE_C_FLAGS       ${CMAKE_C_FLAGS_RELEASE}       )
-        SET(CMAKE_CXX_FLAGS     ${CMAKE_CXX_FLAGS_RELEASE}     )
-        SET(CMAKE_Fortran_FLAGS ${CMAKE_Fortran_FLAGS_RELEASE} )
-    ELSEIF ( ${CMAKE_BUILD_TYPE} STREQUAL "RelWithDebInfo" OR ${CMAKE_BUILD_TYPE} STREQUAL "RELWITHDEBINFO")
-        SET(CMAKE_C_FLAGS       "-g ${CMAKE_C_FLAGS_RELEASE}"       )
-        SET(CMAKE_CXX_FLAGS     "-g ${CMAKE_CXX_FLAGS_RELEASE}"     )
-        SET(CMAKE_Fortran_FLAGS "-g ${CMAKE_Fortran_FLAGS_RELEASE}" )
-    ELSE()
-        MESSAGE(FATAL_ERROR "Unknown value for CMAKE_BUILD_TYPE = ${CMAKE_BUILD_TYPE}")
     ENDIF()
     # Set the behavior of GLIBCXX flags
     CHECK_ENABLE_FLAG( ENABLE_GXX_DEBUG 0 )
@@ -481,6 +477,16 @@ MACRO( SET_COMPILER_FLAGS )
         # Default
         SET( DISABLE_GXX_DEBUG ON )
     ENDIF()
+    # Set debug definitions
+    IF ( ${CMAKE_BUILD_TYPE} STREQUAL "Debug" )
+        SET( CMAKE_C_FLAGS_DEBUG   " ${CMAKE_C_FLAGS_DEBUG}   -DDEBUG -D_DEBUG" )
+        SET( CMAKE_CXX_FLAGS_DEBUG " ${CMAKE_CXX_FLAGS_DEBUG} -DDEBUG -D_DEBUG" )        
+    ENDIF()
+    # Save the debug/release specific flags to the cache
+    SET( CMAKE_C_FLAGS_DEBUG     "${CMAKE_C_FLAGS_DEBUG}"     CACHE STRING "Debug flags"   FORCE)
+    SET( CMAKE_C_FLAGS_RELEASE   "${CMAKE_C_FLAGS_RELEASE}"   CACHE STRING "Release flags" FORCE)
+    SET( CMAKE_CXX_FLAGS_DEBUG   "${CMAKE_CXX_FLAGS_DEBUG}"   CACHE STRING "Debug flags"   FORCE)
+    SET( CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}" CACHE STRING "Release flags" FORCE)
     # Add the user flags
     ADD_USER_FLAGS()
     # Add the c++ standard flags
@@ -492,8 +498,8 @@ MACRO( SET_COMPILER_FLAGS )
     CHECK_CXX_COMPILER_FLAG( "${CMAKE_CXX_FLAGS}" CHECK_CXX_FLAGS )
     IF ( ( NOT CHECK_C_FLAGS ) OR ( NOT CHECK_CXX_FLAGS ) )
         IF ( USING_CRAY )
-            MESSAGE(WARNING "Invalid C/CXX flags detected:"
-                "C flags: ${CMAKE_C_FLAGS}" "CXX flags: ${CMAKE_CXX_FLAGS}" )
+            MESSAGE(WARNING "Invalid C/CXX flags detected:\n"
+                "C flags: ${CMAKE_C_FLAGS}\n" "CXX flags: ${CMAKE_CXX_FLAGS}\n" )
         ENDIF()
     ENDIF()
 ENDMACRO()
@@ -621,6 +627,9 @@ ENDMACRO()
 # Check if we want to keep the test
 FUNCTION( KEEP_TEST RESULT )
     SET( ${RESULT} 1 PARENT_SCOPE )
+    IF ( ONLY_BUILD_DOCS )
+        SET( ${RESULT} 0 PARENT_SCOPE )
+    ENDIF()
 ENDFUNCTION()
 
 
@@ -681,6 +690,7 @@ MACRO( CREATE_TEST_NAME TEST ${ARGN} )
         SET( TESTNAME "${TESTNAME}--${tmp}")
     endforeach()
     # STRING(REGEX REPLACE "--" "-" TESTNAME ${TESTNAME} )
+    SET( LAST_TESTNAME ${TESTNAME} PARENT_SCOPE )
 ENDMACRO()
 
 
@@ -899,21 +909,20 @@ ENDMACRO()
 FUNCTION( ADD_MATLAB_MEX MEXFILE )
     # Set the MEX compiler and default link flags
     IF ( USING_MSVC )
-        #SET(MEX_FLAGS ${MEX_FLAGS} "LINKFLAGS=\"/NODEFAULTLIB:MSVSRT.lib\"" )
-        SET( MEX "\"${MATLAB_DIRECTORY}/sys/perl/win32/bin/perl.exe\" \"${MATLAB_DIRECTORY}/bin/mex.pl\"" )
+        SET( MEX "\"${MATLAB_DIRECTORY}/bin/mex.bat\"" )
         SET( MEX_LDFLAGS ${MEX_LDFLAGS} -L${CMAKE_CURRENT_BINARY_DIR}/.. 
             -L${CMAKE_CURRENT_BINARY_DIR}/../Debug -L${CMAKE_CURRENT_BINARY_DIR} )
-        FOREACH( rpath ${CMAKE_INSTALL_RPATH} )
-            #SET( MEX_LDFLAGS ${MEX_LDFLAGS} "-Wl,-rpath,${rpath},--no-undefined" )
-        ENDFOREACH()
     ELSE()
-        SET( MEX mex )
+        SET( MEX "${MATLAB_DIRECTORY}/bin/mex" )
+        SET( MEX_LDFLAGS ${MEX_LDFLAGS} )
+    ENDIF()
+    SET( MEX_LDFLAGS ${MEX_LDFLAGS} -L${CMAKE_CURRENT_BINARY_DIR}/.. ${SYSTEM_LDFLAGS} )
+    SET( MEX_LIBS    ${MEX_LIBS}    ${SYSTEM_LIBS}    )
+    IF ( NOT USING_MSVC )
         FOREACH( rpath ${CMAKE_INSTALL_RPATH} )
             SET( MEX_LDFLAGS ${MEX_LDFLAGS} "-Wl,-rpath,${rpath},--no-undefined" )
         ENDFOREACH()
     ENDIF()
-    SET( MEX_LDFLAGS ${MEX_LDFLAGS} -L${CMAKE_CURRENT_BINARY_DIR}/.. ${SYSTEM_LDFLAGS} )
-    SET( MEX_LIBS    ${MEX_LIBS}    ${SYSTEM_LIBS}    )
     # Create the mex comamnd
     STRING(REGEX REPLACE "[.]cpp" "" TARGET ${MEXFILE})
     STRING(REGEX REPLACE "[.]c"   "" TARGET ${TARGET})
diff --git a/tests/TestBlobIdentify.cpp b/tests/TestBlobIdentify.cpp
index f79d1c81..3d7b8480 100644
--- a/tests/TestBlobIdentify.cpp
+++ b/tests/TestBlobIdentify.cpp
@@ -24,10 +24,21 @@ inline double rand2()
 
 // Test if all ranks agree on a value
 bool allAgree( int x, MPI_Comm comm ) {
-    int min, max;
-    MPI_Allreduce(&x,&min,1,MPI_INT,MPI_MIN,comm);
-    MPI_Allreduce(&x,&max,1,MPI_INT,MPI_MAX,comm);
-    return min==max;
+    int x2 = x;
+    MPI_Bcast(&x2,1,MPI_INT,0,comm);
+    int diff = x==x2 ? 0:1;
+    int diff2 = 0;
+    MPI_Allreduce(&diff,&diff2,1,MPI_INT,MPI_SUM,comm);
+    return diff2==0;
+}
+template<class T>
+bool allAgree( const std::vector<T>& x, MPI_Comm comm ) {
+    std::vector<T> x2 = x;
+    MPI_Bcast(&x2[0],x.size()*sizeof(T)/sizeof(int),MPI_INT,0,comm);
+    int diff = x==x2 ? 0:1;
+    int diff2 = 0;
+    MPI_Allreduce(&diff,&diff2,1,MPI_INT,MPI_SUM,comm);
+    return diff2==0;
 }
 
 
@@ -247,7 +258,7 @@ int main(int argc, char **argv)
             N_errors++;
         }
         // Identify the blob maps and renumber the ids
-        ID_map_struct map = computeIDMap(GlobalBlobID,GlobalBlobID2,comm);
+        ID_map_struct map = computeIDMap(nx,ny,nz,GlobalBlobID,GlobalBlobID2,comm);
         std::swap(GlobalBlobID,GlobalBlobID2);
         std::vector<BlobIDType> new_list;
         getNewIDs(map,id_max,new_list);
@@ -311,18 +322,28 @@ int main(int argc, char **argv)
         IntArray GlobalBlobID2;
         int nblobs2 = ComputeGlobalBlobIDs(nx,ny,nz,rank_info,Phase,SignDist,vF,vS,GlobalBlobID2,comm);
         // Identify the blob maps and renumber the ids
-        ID_map_struct map = computeIDMap(GlobalBlobID,GlobalBlobID2,comm);
+        ID_map_struct map = computeIDMap(nx,ny,nz,GlobalBlobID,GlobalBlobID2,comm);
         std::swap(GlobalBlobID,GlobalBlobID2);
         std::vector<BlobIDType> new_list;
         getNewIDs(map,id_max,new_list);
-        renumberIDs(new_list,GlobalBlobID);
-        writeIDMap(map,save_it,"lbpm_id_map.txt");
+        // Check id_max
         if ( !allAgree(id_max,comm) ) {
             if ( rank==0 )
                 printf("All ranks do not agree on id_max\n");
             N_errors++;
             break;
         }
+        // Check that the new id list matches on all ranks
+        if ( !allAgree(new_list,comm) ) {
+            if ( rank==0 )
+                printf("All ranks do not agree on new_list\n");
+            N_errors++;
+            break;
+        }
+        // Renumber the ids and write the map
+        renumberIDs(new_list,GlobalBlobID);
+        writeIDMap(map,save_it,"lbpm_id_map.txt");
+        // Check the number of blobs in the map
         int N1 = 0;
         int N2 = 0;
         if ( rank==0 ) {
diff --git a/tests/lbpm_color_simulator.h b/tests/lbpm_color_simulator.h
index 53a3766d..45cd1507 100644
--- a/tests/lbpm_color_simulator.h
+++ b/tests/lbpm_color_simulator.h
@@ -96,7 +96,7 @@ public:
         if ( last_id!=NULL ) {
             // Compute the timestep-timestep map
             const IntArray& old_ids = last_id->second;
-            ID_map_struct map = computeIDMap(old_ids,ids,newcomm);
+            ID_map_struct map = computeIDMap(Nx,Ny,Nz,old_ids,ids,newcomm);
             // Renumber the current timestep's ids
             getNewIDs(map,max_id,*new_list);
             renumberIDs(*new_list,new_id->second);