diff --git a/CMakeLists.txt b/CMakeLists.txt
index f29160e8..5145bb79 100755
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -57,7 +57,7 @@ ADD_CUSTOM_TARGET( doc )
 ADD_CUSTOM_TARGET( latex_docs )
 ADD_CUSTOM_TARGET( build-test )
 ADD_CUSTOM_TARGET( check COMMAND  make test  )
-ADD_DISTCLEAN()
+ADD_DISTCLEAN( tests liblbpm-wia.* cpu gpu example common visit IO )
 
 
 # Check the compile mode and compile flags
@@ -103,6 +103,8 @@ IF ( NOT ONLY_BUILD_DOCS )
     CONFIGURE_LBPM()
     CONFIGURE_TIMER( 0 "${${PROJ}_INSTALL_DIR}/null_timer" )
     CONFIGURE_LINE_COVERAGE()
+    INCLUDE( "${CMAKE_CURRENT_SOURCE_DIR}/cmake/SharedPtr.cmake" )
+    CONFIGURE_SHARED_PTR( "${LBPM_INSTALL_DIR}/include" "" )
 ENDIF()
 
 
diff --git a/IO/Mesh.cpp b/IO/Mesh.cpp
index 3529e5ff..1ca81a28 100644
--- a/IO/Mesh.cpp
+++ b/IO/Mesh.cpp
@@ -2,7 +2,7 @@
 #include "common/Utilities.h"
 
 #include <limits>
-
+#include <stdint.h>
 
 namespace IO {
 
@@ -45,7 +45,7 @@ PointList::~PointList( )
 size_t PointList::numberPointsVar( VariableType type ) const
 {
     size_t N = 0;
-    if ( type == VariableType::NodeVariable )
+    if ( type == NodeVariable )
         N = points.size();
     return N;
 }
@@ -121,9 +121,9 @@ TriList::~TriList( )
 size_t TriList::numberPointsVar( VariableType type ) const
 {
     size_t N = 0;
-    if ( type==VariableType::NodeVariable )
+    if ( type==NodeVariable )
         N = 3*A.size();
-    else if ( type==VariableType::SurfaceVariable ||  type==VariableType::VolumeVariable )
+    else if ( type==SurfaceVariable ||  type==VolumeVariable )
         N = A.size();
     return N;
 }
@@ -192,7 +192,7 @@ TriMesh::TriMesh( size_t N_tri, size_t N_point )
     B.resize(N_tri,-1);
     C.resize(N_tri,-1);
 }
-TriMesh::TriMesh( size_t N_tri, std::shared_ptr<PointList> points )
+TriMesh::TriMesh( size_t N_tri, shared_ptr<PointList> points )
 {
     vertices = points;
     A.resize(N_tri,-1);
@@ -226,9 +226,9 @@ TriMesh::~TriMesh( )
 size_t TriMesh::numberPointsVar( VariableType type ) const
 {
     size_t N = 0;
-    if ( type==VariableType::NodeVariable )
+    if ( type==NodeVariable )
         N = vertices->points.size();
-    else if ( type==VariableType::SurfaceVariable ||  type==VariableType::VolumeVariable )
+    else if ( type==SurfaceVariable ||  type==VolumeVariable )
         N = A.size();
     return N;
 }
@@ -293,45 +293,45 @@ void TriMesh::unpack( const std::pair<size_t,void*>& data_in )
 /****************************************************
 * Converters                                        *
 ****************************************************/
-std::shared_ptr<PointList> getPointList( std::shared_ptr<Mesh> mesh )
+shared_ptr<PointList> getPointList( shared_ptr<Mesh> mesh )
 {
-    return std::dynamic_pointer_cast<PointList>(mesh);
+    return dynamic_pointer_cast<PointList>(mesh);
 }
-std::shared_ptr<TriMesh> getTriMesh( std::shared_ptr<Mesh> mesh )
+shared_ptr<TriMesh> getTriMesh( shared_ptr<Mesh> mesh )
 {
-    std::shared_ptr<TriMesh> mesh2;
-    if ( std::dynamic_pointer_cast<TriMesh>(mesh) != NULL ) {
-        mesh2 = std::dynamic_pointer_cast<TriMesh>(mesh);
-    } else if ( std::dynamic_pointer_cast<TriList>(mesh) != NULL ) {
-        std::shared_ptr<TriList> trilist = std::dynamic_pointer_cast<TriList>(mesh);
-        ASSERT(trilist!=NULL);
+    shared_ptr<TriMesh> mesh2;
+    if ( dynamic_pointer_cast<TriMesh>(mesh).get() != NULL ) {
+        mesh2 = dynamic_pointer_cast<TriMesh>(mesh);
+    } else if ( dynamic_pointer_cast<TriList>(mesh).get() != NULL ) {
+        shared_ptr<TriList> trilist = dynamic_pointer_cast<TriList>(mesh);
+        ASSERT(trilist.get()!=NULL);
         mesh2.reset( new TriMesh(*trilist) );
     }
     return mesh2;
 }
-std::shared_ptr<TriList> getTriList( std::shared_ptr<Mesh> mesh )
+shared_ptr<TriList> getTriList( shared_ptr<Mesh> mesh )
 {
-    std::shared_ptr<TriList> mesh2;
-    if ( std::dynamic_pointer_cast<TriList>(mesh) != NULL ) {
-        mesh2 = std::dynamic_pointer_cast<TriList>(mesh);
-    } else if ( std::dynamic_pointer_cast<TriMesh>(mesh) != NULL ) {
-        std::shared_ptr<TriMesh> trimesh = std::dynamic_pointer_cast<TriMesh>(mesh);
-        ASSERT(trimesh!=NULL);
+    shared_ptr<TriList> mesh2;
+    if ( dynamic_pointer_cast<TriList>(mesh).get() != NULL ) {
+        mesh2 = dynamic_pointer_cast<TriList>(mesh);
+    } else if ( dynamic_pointer_cast<TriMesh>(mesh).get() != NULL ) {
+        shared_ptr<TriMesh> trimesh = dynamic_pointer_cast<TriMesh>(mesh);
+        ASSERT(trimesh.get()!=NULL);
         mesh2.reset( new TriList(*trimesh) );
     }
     return mesh2;
 }
-std::shared_ptr<const PointList> getPointList( std::shared_ptr<const Mesh> mesh )
+shared_ptr<const PointList> getPointList( shared_ptr<const Mesh> mesh )
 {
-    return getPointList( std::const_pointer_cast<Mesh>(mesh) );
+    return getPointList( const_pointer_cast<Mesh>(mesh) );
 }
-std::shared_ptr<const TriMesh> getTriMesh( std::shared_ptr<const Mesh> mesh )
+shared_ptr<const TriMesh> getTriMesh( shared_ptr<const Mesh> mesh )
 {
-    return getTriMesh( std::const_pointer_cast<Mesh>(mesh) );
+    return getTriMesh( const_pointer_cast<Mesh>(mesh) );
 }
-std::shared_ptr<const TriList> getTriList( std::shared_ptr<const Mesh> mesh )
+shared_ptr<const TriList> getTriList( shared_ptr<const Mesh> mesh )
 {
-    return getTriList( std::const_pointer_cast<Mesh>(mesh) );
+    return getTriList( const_pointer_cast<Mesh>(mesh) );
 }
 
 
diff --git a/IO/Mesh.h b/IO/Mesh.h
index 0b3c7ad5..25933ddf 100644
--- a/IO/Mesh.h
+++ b/IO/Mesh.h
@@ -3,17 +3,19 @@
 
 #include <iostream>
 #include <string.h>
-#include <memory>
 #include <vector>
 
 #include "common/PointList.h"
+#include "shared_ptr.h"
+
 
 
 namespace IO {
 
 
 //! Possible variable types
-enum class VariableType : unsigned char { NodeVariable=1, EdgeVariable=2, SurfaceVariable=2, VolumeVariable=3, Null=0 };
+//enum class VariableType : unsigned char { NodeVariable=1, EdgeVariable=2, SurfaceVariable=2, VolumeVariable=3, Null=0 };
+enum VariableType { NodeVariable=1, EdgeVariable=2, SurfaceVariable=2, VolumeVariable=3, NullVariable=0 };
 
 
 /*! \class Mesh
@@ -65,42 +67,10 @@ public:
 };
 
 
-/*! \class TriMesh
-    \brief A class used to hold a list of trianges specified by their vertex number and list of coordiantes
-*/
-class TriList;
-class TriMesh: public Mesh
-{
-public:
-    //! TriMesh constructor
-    TriMesh();
-    //! Constructor for Nt triangles and Np points
-    TriMesh( size_t N_tri, size_t N_point );
-    //! Constructor for Nt triangles and the given points
-    TriMesh( size_t N_tri, std::shared_ptr<PointList> points );
-    //! Constructor from TriList
-    TriMesh( const TriList& );
-    //! Destructor
-    virtual ~TriMesh();
-    //! Mesh class name
-    virtual std::string className() const { return "TriMesh"; }
-    //! Number of points for the given variable type
-    virtual size_t numberPointsVar( VariableType type ) const;
-    //! Pack the data
-    virtual std::pair<size_t,void*> pack( int level ) const;
-    //! Unpack the data
-    virtual void unpack( const std::pair<size_t,void*>& data );
-public:
-    std::shared_ptr<PointList> vertices;    //!< List of verticies
-    std::vector<int>    A;                  //!< First vertex
-    std::vector<int>    B;                  //!< Second vertex
-    std::vector<int>    C;                  //!< Third vertex
-};
-
-
 /*! \class TriList
     \brief A class used to hold a list of triangles specified by their vertex coordinates
 */
+class TriMesh;
 class TriList: public Mesh
 {
 public:
@@ -127,6 +97,38 @@ public:
 };
 
 
+/*! \class TriMesh
+    \brief A class used to hold a list of trianges specified by their vertex number and list of coordiantes
+*/
+class TriMesh: public Mesh
+{
+public:
+    //! TriMesh constructor
+    TriMesh();
+    //! Constructor for Nt triangles and Np points
+    TriMesh( size_t N_tri, size_t N_point );
+    //! Constructor for Nt triangles and the given points
+    TriMesh( size_t N_tri, shared_ptr<PointList> points );
+    //! Constructor from TriList
+    TriMesh( const TriList& );
+    //! Destructor
+    virtual ~TriMesh();
+    //! Mesh class name
+    virtual std::string className() const { return "TriMesh"; }
+    //! Number of points for the given variable type
+    virtual size_t numberPointsVar( VariableType type ) const;
+    //! Pack the data
+    virtual std::pair<size_t,void*> pack( int level ) const;
+    //! Unpack the data
+    virtual void unpack( const std::pair<size_t,void*>& data );
+public:
+    shared_ptr<PointList> vertices;    //!< List of verticies
+    std::vector<int>    A;                  //!< First vertex
+    std::vector<int>    B;                  //!< Second vertex
+    std::vector<int>    C;                  //!< Third vertex
+};
+
+
 
 /*! \class Variable
     \brief A base class fore variables
@@ -140,7 +142,7 @@ public:
     std::string name;           //!< Variable name
     std::vector<double> data;   //!< Variable data
     //! Empty constructor
-    Variable(): type(VariableType::Null) {}
+    Variable(): type(NullVariable) {}
     //! Destructor
     virtual ~Variable() {}
 protected:
@@ -156,18 +158,18 @@ protected:
 */
 struct MeshDataStruct {
     std::string             meshName;
-    std::shared_ptr<Mesh>   mesh;
-    std::vector<std::shared_ptr<Variable> >  vars;
+    shared_ptr<Mesh>   mesh;
+    std::vector<shared_ptr<Variable> >  vars;
 };
 
 
 //! Convert the mesh to a TriMesh (will return NULL if this is invalid)
-std::shared_ptr<PointList> getPointList( std::shared_ptr<Mesh> mesh );
-std::shared_ptr<TriMesh> getTriMesh( std::shared_ptr<Mesh> mesh );
-std::shared_ptr<TriList> getTriList( std::shared_ptr<Mesh> mesh );
-std::shared_ptr<const PointList> getPointList( std::shared_ptr<const Mesh> mesh );
-std::shared_ptr<const TriMesh> getTriMesh( std::shared_ptr<const Mesh> mesh );
-std::shared_ptr<const TriList> getTriList( std::shared_ptr<const Mesh> mesh );
+shared_ptr<PointList> getPointList( shared_ptr<Mesh> mesh );
+shared_ptr<TriMesh> getTriMesh( shared_ptr<Mesh> mesh );
+shared_ptr<TriList> getTriList( shared_ptr<Mesh> mesh );
+shared_ptr<const PointList> getPointList( shared_ptr<const Mesh> mesh );
+shared_ptr<const TriMesh> getTriMesh( shared_ptr<const Mesh> mesh );
+shared_ptr<const TriList> getTriList( shared_ptr<const Mesh> mesh );
 
 
 } // IO namespace
diff --git a/IO/MeshDatabase.cpp b/IO/MeshDatabase.cpp
index 6a30d13f..74d51fd9 100644
--- a/IO/MeshDatabase.cpp
+++ b/IO/MeshDatabase.cpp
@@ -7,6 +7,8 @@
 #include <vector>
 #include <map>
 #include <set>
+#include <cstdio>
+
 #include <ProfilerApp.h>
 
 
@@ -397,13 +399,13 @@ std::vector<MeshDatabase> read( const std::string& filename )
 
 
 // Return the mesh type
-IO::MeshType meshType( std::shared_ptr<IO::Mesh> mesh )
+IO::MeshType meshType( shared_ptr<IO::Mesh> mesh )
 {
-    IO::MeshType type = IO::MeshType::Unknown;
-    if ( std::dynamic_pointer_cast<IO::PointList>(mesh)!=NULL ) {
-        type = IO::MeshType::PointMesh;
-    } else if ( std::dynamic_pointer_cast<IO::TriList>(mesh)!=NULL || std::dynamic_pointer_cast<IO::TriMesh>(mesh)!=NULL ) {
-        type = IO::MeshType::SurfaceMesh;
+    IO::MeshType type = IO::Unknown;
+    if ( dynamic_pointer_cast<IO::PointList>(mesh).get()!=NULL ) {
+        type = IO::PointMesh;
+    } else if ( dynamic_pointer_cast<IO::TriList>(mesh).get()!=NULL || dynamic_pointer_cast<IO::TriMesh>(mesh).get()!=NULL ) {
+        type = IO::SurfaceMesh;
     } else {
         ERROR("Unknown mesh");
     }
diff --git a/IO/MeshDatabase.h b/IO/MeshDatabase.h
index 5587be9d..abf74f23 100644
--- a/IO/MeshDatabase.h
+++ b/IO/MeshDatabase.h
@@ -3,10 +3,10 @@
 
 #include "IO/Mesh.h" 
 #include "common/MPI_Helpers.h"
+#include "shared_ptr.h"
 
 #include <iostream>
 #include <string.h>
-#include <memory>
 #include <vector>
 #include <map>
 
@@ -17,7 +17,8 @@ class Mesh;
 
 
 //! Enum to identify mesh type
-enum class MeshType : char { PointMesh=1, SurfaceMesh=2, VolumeMesh=3, Unknown=-1 };
+//enum class MeshType : char { PointMesh=1, SurfaceMesh=2, VolumeMesh=3, Unknown=-1 };
+enum MeshType { PointMesh=1, SurfaceMesh=2, VolumeMesh=3, Unknown=-1 };
 
 
 //! Helper struct for containing offsets for the mesh info
@@ -80,7 +81,7 @@ std::vector<MeshDatabase> read( const std::string& filename );
 
 
 //! Return the mesh type
-IO::MeshType meshType( std::shared_ptr<IO::Mesh> mesh );
+IO::MeshType meshType( shared_ptr<IO::Mesh> mesh );
 
 
 } // IO namespace
diff --git a/IO/PIO.cpp b/IO/PIO.cpp
index 34cabeca..6c6ece2d 100644
--- a/IO/PIO.cpp
+++ b/IO/PIO.cpp
@@ -135,10 +135,11 @@ void ParallelStreamBuffer::reserve( size_t size )
         }
     }
 }
-std::streamsize ParallelStreamBuffer::xsputn( const std::string &text, std::streamsize n )
+std::streamsize ParallelStreamBuffer::xsputn( const char* text, std::streamsize n )
 {
     reserve(d_size+n);
-    memcpy(&d_buffer[d_size],text.c_str(),text.size());
+    memcpy(&d_buffer[d_size],text,n);
+    d_size += n;
     if ( text[n-1]==0 || text[n-1]==10 ) { sync(); }
     return n;
 }
diff --git a/IO/PIO.h b/IO/PIO.h
index e21d7624..b6d8b103 100644
--- a/IO/PIO.h
+++ b/IO/PIO.h
@@ -73,7 +73,7 @@ public:
      * Write the specified number of characters into the output stream (called
      * from streambuf).
     */
-    virtual std::streamsize xsputn(const std::string &text, std::streamsize n);
+    virtual std::streamsize xsputn(const char* text, std::streamsize n);
 
     /*!
      * Write an overflow character into the parallel buffer (called from
diff --git a/IO/Reader.cpp b/IO/Reader.cpp
index 9f3133bc..d6144a87 100644
--- a/IO/Reader.cpp
+++ b/IO/Reader.cpp
@@ -7,11 +7,18 @@
 #include <ProfilerApp.h>
 #include <iostream>
 #include <string.h>
-#include <memory>
 #include <vector>
 #include <map>
+#include <cstdio>
 
 
+// Inline function to read line without a return argument
+static inline void fgetl( char * str, int num, FILE * stream )
+{
+    char* ptr = fgets( str, num, stream );
+    if ( 0 ) {char *temp = (char *)&ptr; temp++;}
+}
+
 
 // List the timesteps in the given directors (dumps.LBPM)
 std::vector<std::string> IO::readTimesteps( const std::string& filename )
@@ -44,11 +51,11 @@ std::vector<IO::MeshDatabase> IO::getMeshList( const std::string& path, const st
 
 
 // Read the given mesh domain
-std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::string& timestep, 
+shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::string& timestep, 
     const IO::MeshDatabase& meshDatabase, int domain )
 {
     PROFILE_START("getMesh");
-    std::shared_ptr<IO::Mesh> mesh;
+    shared_ptr<IO::Mesh> mesh;
     if ( meshDatabase.format==1 ) {
         // Old format (binary doubles)
         std::string filename = path + "/" + timestep + "/" + meshDatabase.domains[domain].file;
@@ -63,9 +70,9 @@ std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::strin
         fclose(fid);
         if ( count%3 != 0 )
             ERROR("Error reading file");
-        if ( meshDatabase.type==IO::MeshType::PointMesh ) {
+        if ( meshDatabase.type==IO::PointMesh ) {
             size_t N = count/3;
-            std::shared_ptr<PointList> pointlist( new PointList(N) );
+            shared_ptr<PointList> pointlist( new PointList(N) );
             std::vector<Point>& P = pointlist->points;
             for (size_t i=0; i<N; i++) {
                 P[i].x = data[3*i+0];
@@ -73,11 +80,11 @@ std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::strin
                 P[i].z = data[3*i+2];
             }
             mesh = pointlist;
-        } else if ( meshDatabase.type==IO::MeshType::SurfaceMesh ) {
+        } else if ( meshDatabase.type==IO::SurfaceMesh ) {
             if ( count%9 != 0 )
                 ERROR("Error reading file (2)");
             size_t N_tri = count/9;
-            std::shared_ptr<TriList> trilist( new TriList(N_tri) );
+            shared_ptr<TriList> trilist( new TriList(N_tri) );
             std::vector<Point>& A = trilist->A;
             std::vector<Point>& B = trilist->B;
             std::vector<Point>& C = trilist->C;
@@ -103,7 +110,7 @@ std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::strin
         FILE *fid = fopen(filename.c_str(),"rb");
         fseek(fid,database.offset,SEEK_SET);
         char line[1000];
-        std::fgets(line,1000,fid);
+        fgetl(line,1000,fid);
         size_t i1 = find(line,':');
         size_t i2 = find(&line[i1+1],':')+i1+1;
         size_t bytes = atol(&line[i2+1]);
@@ -131,20 +138,20 @@ std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::strin
 
 
 // Read the given variable for the given mesh domain
-std::shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const std::string& timestep, 
+shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const std::string& timestep, 
     const MeshDatabase& meshDatabase, int domain, const std::string& variable )
 {
     std::pair<std::string,std::string> key(meshDatabase.domains[domain].name,variable);
     std::map<std::pair<std::string,std::string>,DatabaseEntry>::const_iterator it;
     it = meshDatabase.variable_data.find(key);
     if ( it==meshDatabase.variable_data.end() )
-        return std::shared_ptr<IO::Variable>();
+        return shared_ptr<IO::Variable>();
     const DatabaseEntry& database = it->second;
     std::string filename = path + "/" + timestep + "/" + database.file;
     FILE *fid = fopen(filename.c_str(),"rb");
     fseek(fid,database.offset,SEEK_SET);
     char line[1000];
-    std::fgets(line,1000,fid);
+    fgetl(line,1000,fid);
     size_t i1 = find(line,':');
     size_t i2 = find(&line[i1+1],':')+i1+1;
     std::vector<std::string> values = splitList(&line[i2+1],',');
@@ -158,12 +165,14 @@ std::shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const st
     size_t count = fread(data,1,bytes,fid);
     fclose(fid);
     ASSERT(count==bytes);
-    std::shared_ptr<IO::Variable> var( new IO::Variable() );
+    shared_ptr<IO::Variable> var( new IO::Variable() );
     var->dim = dim;
     var->type = static_cast<IO::VariableType>(type);
     var->name = variable;
     var->data.resize(N);
-    double *var_data = var->data.data();
+    double *var_data = NULL;
+    if ( !var->data.empty() )
+        var_data = &var->data[0];
     if ( precision=="double" ) {
         memcpy(var_data,data,bytes);
     } else {
diff --git a/IO/Reader.h b/IO/Reader.h
index 1a8adead..510cdfdb 100644
--- a/IO/Reader.h
+++ b/IO/Reader.h
@@ -3,11 +3,11 @@
 
 #include <iostream>
 #include <string.h>
-#include <memory>
 #include <vector>
 
 #include "IO/Mesh.h"
 #include "IO/MeshDatabase.h"
+#include "shared_ptr.h"
 
 
 namespace IO {
@@ -22,12 +22,12 @@ std::vector<IO::MeshDatabase> getMeshList( const std::string& path, const std::s
 
 
 //! Read the given mesh domain
-std::shared_ptr<IO::Mesh> getMesh( const std::string& path, const std::string& timestep, 
+shared_ptr<IO::Mesh> getMesh( const std::string& path, const std::string& timestep, 
     const MeshDatabase& meshDatabase, int domain );
 
 
 //! Read the given mesh domain
-std::shared_ptr<IO::Variable> getVariable( const std::string& path, const std::string& timestep, 
+shared_ptr<IO::Variable> getVariable( const std::string& path, const std::string& timestep, 
     const MeshDatabase& meshDatabase, int domain, const std::string& variable );
 
 
diff --git a/IO/Writer.cpp b/IO/Writer.cpp
index 0d0c74d3..455f13a6 100644
--- a/IO/Writer.cpp
+++ b/IO/Writer.cpp
@@ -3,6 +3,7 @@
 #include "IO/IOHelpers.h"
 #include "common/MPI_Helpers.h"
 #include "common/Utilities.h"
+#include "shared_ptr.h"
 
 #include <sys/stat.h>
 #include <algorithm>
@@ -25,7 +26,7 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
 	    sprintf(fullpath,"%s/%s",path,filename);
         FILE *fid = fopen(fullpath,"wb");
         INSIST(fid!=NULL,std::string("Error opening file: ")+fullpath);
-        std::shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
+        shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
         IO::MeshDatabase mesh_entry;
         mesh_entry.name = meshData[i].meshName;
         mesh_entry.type = meshType(mesh);
@@ -41,18 +42,18 @@ static std::vector<IO::MeshDatabase> writeMeshesOrigFormat( const std::vector<IO
             //for (size_t j=0; j<meshData[i].vars.size(); j++)
             //    mesh_entry.variables.push_back( meshData[i].vars[j]->name );
         }
-        if ( std::dynamic_pointer_cast<IO::PointList>(mesh)!=NULL ) {
+        if ( dynamic_pointer_cast<IO::PointList>(mesh).get()!=NULL ) {
             // List of points
-            std::shared_ptr<IO::PointList> pointlist = std::dynamic_pointer_cast<IO::PointList>(mesh);
+            shared_ptr<IO::PointList> pointlist = dynamic_pointer_cast<IO::PointList>(mesh);
             const std::vector<Point>& P = pointlist->points;
             for (size_t i=0; i<P.size(); i++) {
                 double x[3];
                 x[0] = P[i].x;  x[1] = P[i].y;  x[2] = P[i].z;
                 fwrite(x,sizeof(double),3,fid);
             }
-        } else if ( std::dynamic_pointer_cast<IO::TriList>(mesh)!=NULL || std::dynamic_pointer_cast<IO::TriMesh>(mesh)!=NULL ) {
+        } else if ( dynamic_pointer_cast<IO::TriList>(mesh).get()!=NULL || dynamic_pointer_cast<IO::TriMesh>(mesh).get()!=NULL ) {
             // Triangle mesh
-            std::shared_ptr<IO::TriList> trilist = IO::getTriList(mesh);
+            shared_ptr<IO::TriList> trilist = IO::getTriList(mesh);
             const std::vector<Point>& A = trilist->A;
             const std::vector<Point>& B = trilist->B;
             const std::vector<Point>& C = trilist->C;
@@ -117,7 +118,7 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string& filename,
         int type = static_cast<int>(mesh.vars[i]->type);
         size_t N = mesh.vars[i]->data.size();
         const void* data = N==0 ? 0:&mesh.vars[i]->data[0];
-        if ( type == static_cast<int>(IO::VariableType::Null) ) {
+        if ( type == static_cast<int>(IO::NullVariable) ) {
             ERROR("Variable type not set");
         }
         size_t N_mesh = mesh.mesh->numberPointsVar(mesh.vars[i]->type);
@@ -143,7 +144,7 @@ static std::vector<IO::MeshDatabase> writeMeshesNewFormat(
     sprintf(fullpath,"%s/%s",path,filename);
     FILE *fid = fopen(fullpath,"wb");
     for (size_t i=0; i<meshData.size(); i++) {
-        std::shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
+        shared_ptr<IO::Mesh> mesh = meshData[i].mesh;
         meshes_written.push_back( write_domain(fid,filename,meshData[i],format) );
     }
     fclose(fid);
diff --git a/IO/Writer.h b/IO/Writer.h
index 5abb6b2d..2cefdac3 100644
--- a/IO/Writer.h
+++ b/IO/Writer.h
@@ -3,7 +3,6 @@
 
 #include <iostream>
 #include <string.h>
-#include <memory>
 #include <vector>
 
 #include "IO/Mesh.h"
diff --git a/ValgrindSuppresionFile b/ValgrindSuppresionFile
index 222bcc11..34bb2527 100644
--- a/ValgrindSuppresionFile
+++ b/ValgrindSuppresionFile
@@ -1,4 +1,4 @@
-# ACML supressions
+# ACML suppressions
 {
    IdentifyCPUCond
    Memcheck:Cond
@@ -16,6 +16,37 @@
 }
 
 
+# MPI suppressions
+{
+   HYD_pmci_wait_for_completion
+   Memcheck:Leak
+   ...
+   fun:HYD_pmci_wait_for_completion
+   fun:main
+}
+{
+   HYDT_dmxu_poll_wait_for_event
+   Memcheck:Leak
+   ...
+   fun:HYDT_dmxu_poll_wait_for_event
+   fun:main
+}
+{
+   PMPI_Init
+   Memcheck:Leak
+   ...
+   fun:PMPI_Init
+   fun:main
+}
 
 
+# System suppressions
+{
+   expand_dynamic_string_token
+   Memcheck:Cond
+   fun:index
+   fun:expand_dynamic_string_token
+   ...
+}
+
 
diff --git a/cmake/CompareOutput.cmake b/cmake/CompareOutput.cmake
new file mode 100644
index 00000000..85f13006
--- /dev/null
+++ b/cmake/CompareOutput.cmake
@@ -0,0 +1,50 @@
+# This script compares the output of TEST against GOLD, 
+#    ensuring that all lines within GOLD are int TEST.
+# Note that TEST may have additional lines that are not checked
+CMAKE_POLICY(SET CMP0007 OLD)
+
+FILE(READ "${TEST}" output )
+FILE(READ "${GOLD}" sol )
+
+macro(LIST_REPLACE LIST INDEX NEWVALUE)
+    list(INSERT ${LIST} ${INDEX} ${NEWVALUE})
+    MATH(EXPR __INDEX "${INDEX} + 1")
+    list(REMOVE_AT ${LIST} ${__INDEX})
+endmacro(LIST_REPLACE)
+
+# Convert file contents into a CMake list (where each element in the list is one line of the file)
+STRING(REGEX REPLACE ";" "\\\\;" data "${output}")
+STRING(REGEX REPLACE ";" "\\\\;" sol  "${sol}")
+STRING(REGEX REPLACE "\n" ";" data "${data}")
+STRING(REGEX REPLACE "\n" ";" sol "${sol}")
+LIST( LENGTH data N_data )
+LIST( LENGTH sol N_sol )
+MATH( EXPR N_data "${N_data}-1" )
+MATH( EXPR N_sol "${N_sol}-1" )
+FOREACH( index RANGE ${N_data} )
+    LIST(GET data ${index} tmp )
+    STRING(REGEX REPLACE "(\n|\r)" "" tmp "${tmp}")
+    STRING(STRIP "${tmp}" tmp )
+    LIST_REPLACE( data ${index} "${tmp}")
+ENDFOREACH()
+FOREACH( index RANGE ${N_sol} )
+    LIST( GET sol ${index} tmp )
+    STRING(REGEX REPLACE "(\n|\r)" "" tmp "${tmp}")
+    STRING(STRIP "${tmp}" tmp )
+    LIST_REPLACE( sol ${index} "${tmp}")
+ENDFOREACH()
+
+# Check that each line of sol is present in data (and delete it)
+FOREACH( tmp ${sol} )
+    LIST(FIND data "${tmp}" result )
+    IF ( ${result} EQUAL -1 )
+        MESSAGE("Test output:\n${output}\n\n")
+        MESSAGE(FATAL_ERROR "Did not find '${tmp}' in test output\n" )
+    ELSE()
+        LIST(REMOVE_AT data ${result} )
+    ENDIF()
+ENDFOREACH()
+
+# Finished
+MESSAGE( "All lines in ${GOLD} were found in ${TEST}")
+
diff --git a/cmake/SharedPtr.cmake b/cmake/SharedPtr.cmake
new file mode 100644
index 00000000..016128a9
--- /dev/null
+++ b/cmake/SharedPtr.cmake
@@ -0,0 +1,165 @@
+# Create a shared_ptr.h file in the include directory that contains 
+#    a shared_ptr class (hopefully typedef to a compiler basic)
+# Arguements:
+#    INSTALL_DIR - Directory to install shared_ptr.h
+#    NAMESPACE - Namespace to contain the shared_ptr class (may be empty)
+INCLUDE( CheckCXXSourceCompiles )
+FUNCTION( CONFIGURE_SHARED_PTR INSTALL_DIR NAMESPACE )
+    SET( CMAKE_REQUIRED_FLAGS ${CMAKE_CXX_FLAGS} )
+    CHECK_CXX_SOURCE_COMPILES(
+	    "   #include <memory>
+            namespace ${NAMESPACE} { using std::shared_ptr; }
+	        int main() {
+	            ${NAMESPACE}::shared_ptr<int> ptr;
+	            return 0;
+	        }
+	    "
+	    MEMORY_SHARED_PTR )
+    CHECK_CXX_SOURCE_COMPILES(
+	    "   #include <memory>
+            namespace ${NAMESPACE} { using std::tr1::shared_ptr; }
+	        int main() {
+	            ${NAMESPACE}::shared_ptr<int> ptr;
+	            return 0;
+	        }
+	    "
+	    MEMORY_TR1_SHARED_PTR )
+    CHECK_CXX_SOURCE_COMPILES(
+	    "   #include <tr1/memory>
+            namespace  ${NAMESPACE} { using std::tr1::shared_ptr; }
+	        int main() {
+	            ${NAMESPACE}::shared_ptr<int> ptr;
+	            return 0;
+	        }
+	    "
+	    TR1_MEMORY_TR1_SHARED_PTR )
+    GET_DIRECTORY_PROPERTY( dirs INCLUDE_DIRECTORIES )
+    SET( CMAKE_REQUIRED_FLAGS "${CMAKE_CXX_FLAGS}" )
+    SET( CMAKE_REQUIRED_INCLUDES ${dirs} "${BOOST_INCLUDE}" )
+    CHECK_CXX_SOURCE_COMPILES(
+	    "   #include \"boost/shared_ptr.hpp\"
+            namespace  ${NAMESPACE} { using boost::shared_ptr; }
+	        int main() {
+	            ${NAMESPACE}::shared_ptr<int> ptr;
+	            return 0;
+	        }
+	    "
+	    BOOST_SHARED_PTR )
+    WRITE_DUMMY_SHARED_PTR( "${NAMESPACE}" "${CMAKE_CURRENT_BINARY_DIR}/tmp/dummy_shared_ptr.h" )
+    CHECK_CXX_SOURCE_COMPILES(
+	    "   #include <iostream>
+	        #include \"${CMAKE_CURRENT_BINARY_DIR}/tmp/dummy_shared_ptr.h\"
+	        int main() {
+	            ${NAMESPACE}::shared_ptr<int> ptr;
+	            return 0;
+	        }
+	    "
+	    DUMMY_SHARED_PTR )
+    IF ( NOT NAMESPACE )
+        SET( NAMESPACE " " )
+    ENDIF()
+    IF ( BOOST_SHARED_PTR )
+        FILE(WRITE  "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "#include \"boost/shared_ptr.hpp\"\n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "#include \"boost/weak_ptr.hpp\"\n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "#include \"boost/enable_shared_from_this.hpp\"\n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "namespace ${NAMESPACE} {\n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using boost::shared_ptr; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using boost::dynamic_pointer_cast; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using boost::const_pointer_cast; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using boost::weak_ptr; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using boost::enable_shared_from_this; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "}\n")
+    ELSEIF ( MEMORY_SHARED_PTR )
+        IF ( ${NAMESPACE} STREQUAL "std" )
+            FILE(WRITE  "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "#include <memory>\n")
+        ELSE()
+            FILE(WRITE  "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "#include <memory>\n")
+            FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "namespace ${NAMESPACE} {\n")
+            FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::shared_ptr; \n")
+            FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::dynamic_pointer_cast; \n")
+            FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::const_pointer_cast; \n")
+            FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::weak_ptr; \n")
+            FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::enable_shared_from_this; \n")
+            FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "}\n")
+        ENDIF()
+    ELSEIF ( MEMORY_TR1_SHARED_PTR )
+        FILE(WRITE  "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "#include <memory>\n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "namespace ${NAMESPACE} {\n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::tr1::shared_ptr; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::tr1::dynamic_pointer_cast; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::tr1::const_pointer_cast; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::tr1::weak_ptr; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::tr1::enable_shared_from_this; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "}\n")
+    ELSEIF ( TR1_MEMORY_TR1_SHARED_PTR )
+        FILE(WRITE  "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "#include <tr1/memory>\n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "namespace ${NAMESPACE} {\n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::tr1::shared_ptr; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::tr1::dynamic_pointer_cast; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::tr1::const_pointer_cast; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::tr1::weak_ptr; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "    using std::tr1::enable_shared_from_this; \n")
+        FILE(APPEND "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "}\n")
+    ELSEIF ( DUMMY_SHARED_PTR ) 
+        MESSAGE("Warning: No valid shared_ptr found, using dummy shared_ptr" )
+        WRITE_DUMMY_SHARED_PTR( "${NAMESPACE}" "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" )
+    ELSE()
+        MESSAGE(FATAL_ERROR "No shared_ptr availible")
+    ENDIF()
+    EXECUTE_PROCESS( COMMAND ${CMAKE_COMMAND} -E copy_if_different 
+        "${CMAKE_CURRENT_BINARY_DIR}/tmp/shared_ptr.h" "${INSTALL_DIR}/shared_ptr.h" )
+ENDFUNCTION()
+
+
+FUNCTION( WRITE_DUMMY_SHARED_PTR NAMESPACE FILENAME )
+    FILE(WRITE  "${FILENAME}" "#ifndef DUMMY_SHARED_PTR_INC\n")
+    FILE(APPEND "${FILENAME}" "#define DUMMY_SHARED_PTR_INC\n")
+    FILE(APPEND "${FILENAME}" "namespace dummy {\n\n")
+    FILE(APPEND "${FILENAME}" "template<class T> class shared_ptr {\n")
+    FILE(APPEND "${FILENAME}" "public:\n")
+    FILE(APPEND "${FILENAME}" "    shared_ptr( ): obj(NULL), count(NULL) {}\n")
+    FILE(APPEND "${FILENAME}" "    shared_ptr( T *ptr ): obj(ptr), count(NULL) { if (ptr) { count = new int; (*count)=1; } } \n")
+    FILE(APPEND "${FILENAME}" "    shared_ptr( const shared_ptr<T>& rhs ):  \n")
+    FILE(APPEND "${FILENAME}" "        obj(rhs.get()), count(rhs.count) { if ( count!=NULL ) { ++(*count); } } \n")
+    FILE(APPEND "${FILENAME}" "    template<class U> shared_ptr( const shared_ptr<U>& rhs ):  \n")
+    FILE(APPEND "${FILENAME}" "        obj(rhs.get()), count(rhs.count) { if ( count!=NULL ) { ++(*count); } } \n")
+    FILE(APPEND "${FILENAME}" "    shared_ptr& operator=( const shared_ptr<T>& rhs ) { obj=rhs.obj; count=rhs.count; ++(*count); return *this; } \n")
+    FILE(APPEND "${FILENAME}" "    ~shared_ptr( ) { reset(); }\n")
+    FILE(APPEND "${FILENAME}" "    void reset( T *ptr ) { reset(); obj=ptr; count=new int; (*count)=1; }\n")
+    FILE(APPEND "${FILENAME}" "    void reset( void ) { \n")
+    FILE(APPEND "${FILENAME}" "        if ( count!=NULL) { int tmp=--(*count); if ( tmp==0 ) { delete obj; delete count; } } \n")
+    FILE(APPEND "${FILENAME}" "        obj=NULL; count=NULL; \n")
+    FILE(APPEND "${FILENAME}" "    }\n")
+    FILE(APPEND "${FILENAME}" "    T* get( ) const { return obj; } \n")
+    FILE(APPEND "${FILENAME}" "    T* operator->( ) const { return obj; } \n")
+    FILE(APPEND "${FILENAME}" "    const T& operator*( ) const { return *obj; } \n")
+    FILE(APPEND "${FILENAME}" "    bool operator==( const T * rhs ) const { return obj==rhs; } \n")
+    FILE(APPEND "${FILENAME}" "    bool operator!=( const T * rhs ) const { return obj!=rhs; } \n")
+    FILE(APPEND "${FILENAME}" "protected:\n")
+    FILE(APPEND "${FILENAME}" "    T *obj;\n")
+    FILE(APPEND "${FILENAME}" "    volatile int *count;\n")
+    FILE(APPEND "${FILENAME}" "template<class T1, class U> friend shared_ptr<T1> dynamic_pointer_cast( shared_ptr<U> const & );\n")
+    FILE(APPEND "${FILENAME}" "template<class T1, class U> friend shared_ptr<T1> const_pointer_cast( shared_ptr<U> const & );\n")
+    FILE(APPEND "${FILENAME}" "template<class Y> friend class shared_ptr;\n")
+    FILE(APPEND "${FILENAME}" "};\n\n")
+    FILE(APPEND "${FILENAME}" "template<class T, class U> shared_ptr<T> dynamic_pointer_cast( shared_ptr<U> const & rhs ) {\n")
+    FILE(APPEND "${FILENAME}" "    T* obj = dynamic_cast<T*>(rhs.obj);\n")
+    FILE(APPEND "${FILENAME}" "    shared_ptr<T> ptr;\n")
+    FILE(APPEND "${FILENAME}" "    if ( obj!=NULL ) { ptr.obj = obj; ptr.count=rhs.count; ++(*ptr.count); }\n")
+    FILE(APPEND "${FILENAME}" "    return ptr;\n}\n")
+    FILE(APPEND "${FILENAME}" "template<class T, class U> shared_ptr<T> const_pointer_cast( shared_ptr<U> const & rhs ) {\n")
+    FILE(APPEND "${FILENAME}" "    T* obj = const_cast<T*>(rhs.obj);\n")
+    FILE(APPEND "${FILENAME}" "    shared_ptr<T> ptr;\n")
+    FILE(APPEND "${FILENAME}" "    if ( obj!=NULL ) { ptr.obj = obj; ptr.count=rhs.count; ++(*ptr.count); }\n")
+    FILE(APPEND "${FILENAME}" "    return ptr;\n}\n")
+    FILE(APPEND "${FILENAME}" "\n} // namespace dummy\n")
+    FILE(APPEND "${FILENAME}" "\n\n")
+    FILE(APPEND "${FILENAME}" "namespace ${NAMESPACE} {\n")
+    FILE(APPEND "${FILENAME}" "    using dummy::shared_ptr; \n")
+    FILE(APPEND "${FILENAME}" "    using dummy::dynamic_pointer_cast; \n")
+    FILE(APPEND "${FILENAME}" "    using dummy::const_pointer_cast; \n")
+    FILE(APPEND "${FILENAME}" "}\n\n")
+    FILE(APPEND "${FILENAME}" "#endif\n")
+ENDFUNCTION()
+
+
diff --git a/cmake/ctest_script.cmake b/cmake/ctest_script.cmake
index 2cd0b1a6..f654a245 100644
--- a/cmake/ctest_script.cmake
+++ b/cmake/ctest_script.cmake
@@ -6,8 +6,12 @@
 
 # Set platform specific variables
 SITE_NAME( HOSTNAME )
+STRING(REGEX REPLACE "-ext."   "" HOSTNAME "${HOSTNAME}")
+STRING(REGEX REPLACE "-login." "" HOSTNAME "${HOSTNAME}")
 SET( CC                 $ENV{CC}                )
 SET( CXX                $ENV{CXX}               )
+SET( CFLAGS             $ENV{CFLAGS}            )
+SET( CXXFLAGS           $ENV{CXXFLAGS}          )
 SET( MPIEXEC            $ENV{MPIEXEC}           )
 SET( USE_TIMER          "$ENV{USE_TIMER}"       )
 SET( TIMER_DIRECTORY    "$ENV{TIMER_DIRECTORY}" )
@@ -35,6 +39,11 @@ SET( MPI_LINK_FLAGS     $ENV{MPI_LINK_FLAGS}    )
 SET( MPI_LIBRARIES      $ENV{MPI_LIBRARIES}     )
 SET( MPIEXEC            $ENV{MPIEXEC}           )
 SET( BUILD_SERIAL       $ENV{BUILD_SERIAL}      )
+SET( CUDA_FLAGS         $ENV{CUDA_FLAGS}        )
+SET( CUDA_HOST_COMPILER $ENV{CUDA_HOST_COMPILER} )
+SET( SKIP_TESTS         $ENV{SKIP_TESTS}        )
+SET( BUILDNAME_POSTFIX  "$ENV{BUILDNAME_POSTFIX}" )
+SET( LIB_TYPE           "$ENV{LIB_TYPE}"        )
 
 
 # Get the source directory based on the current directory
@@ -47,53 +56,53 @@ ENDIF()
 
 
 # Check that we specified the build type to run
+SET( USE_VALGRIND FALSE )
+SET( RUN_WEEKLY FALSE )
+SET( USE_CUDA FALSE )
+SET( ENABLE_GCOV "false" )
+SET( CTEST_COVERAGE_COMMAND ${COVERAGE_COMMAND} )
 IF( NOT CTEST_SCRIPT_ARG )
     MESSAGE(FATAL_ERROR "No build specified: ctest -S /path/to/script,build (debug/optimized/valgrind")
 ELSEIF( ${CTEST_SCRIPT_ARG} STREQUAL "debug" )
     SET( CTEST_BUILD_NAME "LBPM-WIA-debug" )
     SET( CMAKE_BUILD_TYPE "Debug" )
-    SET( CTEST_COVERAGE_COMMAND ${COVERAGE_COMMAND} )
     SET( ENABLE_GCOV "true" )
-    SET( USE_VALGRIND FALSE )
-    SET( USE_CUDA FALSE )
 ELSEIF( ${CTEST_SCRIPT_ARG} STREQUAL "debug-cuda" )
     SET( CTEST_BUILD_NAME "LBPM-WIA-debug-cuda" )
     SET( CMAKE_BUILD_TYPE "Debug" )
-    SET( CTEST_COVERAGE_COMMAND ${COVERAGE_COMMAND} )
     SET( ENABLE_GCOV "true" )
-    SET( USE_VALGRIND FALSE )
     SET( USE_CUDA TRUE )
 ELSEIF( (${CTEST_SCRIPT_ARG} STREQUAL "optimized") OR (${CTEST_SCRIPT_ARG} STREQUAL "opt") )
     SET( CTEST_BUILD_NAME "LBPM-WIA-opt" )
     SET( CMAKE_BUILD_TYPE "Release" )
-    SET( CTEST_COVERAGE_COMMAND )
-    SET( ENABLE_GCOV "false" )
-    SET( USE_VALGRIND FALSE )
-    SET( USE_CUDA FALSE )
 ELSEIF( (${CTEST_SCRIPT_ARG} STREQUAL "optimized-cuda") OR (${CTEST_SCRIPT_ARG} STREQUAL "opt-cuda") )
     SET( CTEST_BUILD_NAME "LBPM-WIA-opt-cuda" )
     SET( CMAKE_BUILD_TYPE "Release" )
-    SET( CTEST_COVERAGE_COMMAND )
-    SET( ENABLE_GCOV "false" )
-    SET( USE_VALGRIND FALSE )
+    SET( USE_CUDA TRUE )
+ELSEIF( (${CTEST_SCRIPT_ARG} STREQUAL "weekly") )
+    SET( CTEST_BUILD_NAME "LBPM-WIA-weekly" )
+    SET( CMAKE_BUILD_TYPE "Release" )
+    SET( RUN_WEEKLY TRUE )
+ELSEIF( (${CTEST_SCRIPT_ARG} STREQUAL "weekly-cuda") )
+    SET( CTEST_BUILD_NAME "LBPM-WIA-weekly-cuda" )
+    SET( CMAKE_BUILD_TYPE "Release" )
+    SET( RUN_WEEKLY TRUE )
     SET( USE_CUDA TRUE )
 ELSEIF( ${CTEST_SCRIPT_ARG} STREQUAL "valgrind" )
     SET( CTEST_BUILD_NAME "LBPM-WIA-valgrind" )
     SET( CMAKE_BUILD_TYPE "Debug" )
-    SET( CTEST_COVERAGE_COMMAND )
-    SET( ENABLE_GCOV "false" )
     SET( USE_VALGRIND TRUE )
-    SET( USE_CUDA FALSE )
 ELSEIF( ${CTEST_SCRIPT_ARG} STREQUAL "valgrind-cuda" )
     SET( CTEST_BUILD_NAME "LBPM-WIA-valgrind-cuda" )
     SET( CMAKE_BUILD_TYPE "Debug" )
-    SET( CTEST_COVERAGE_COMMAND )
-    SET( ENABLE_GCOV "false" )
     SET( USE_VALGRIND TRUE )
     SET( USE_CUDA TRUE )
 ELSE()
     MESSAGE(FATAL_ERROR "Invalid build (${CTEST_SCRIPT_ARG}): ctest -S /path/to/script,build (debug/opt/valgrind")
 ENDIF()
+IF ( BUILDNAME_POSTFIX )
+    SET( CTEST_BUILD_NAME "${CTEST_BUILD_NAME}-${BUILDNAME_POSTFIX}" )
+ENDIF()
 IF ( NOT CTEST_COVERAGE_COMMAND )
     SET( ENABLE_GCOV "false" )
 ENDIF()
@@ -113,10 +122,10 @@ IF( NOT DEFINED N_PROCS )
     ENDIF()
     # Mac:
     IF(APPLE)
-        find_program(cmd_sys_pro "system_profiler")
+        find_program(cmd_sys_pro "sysctl")
         if(cmd_sys_pro)
-            execute_process(COMMAND ${cmd_sys_pro} OUTPUT_VARIABLE info)
-            STRING(REGEX REPLACE "^.*Total Number of Cores: ([0-9]+).*$" "\\1" N_PROCS "${info}")
+            execute_process(COMMAND ${cmd_sys_pro} hw.physicalcpu OUTPUT_VARIABLE info)
+            STRING(REGEX REPLACE "^.*hw.physicalcpu: ([0-9]+).*$" "\\1" N_PROCS "${info}")
         ENDIF()
     ENDIF()
     # Windows:
@@ -143,6 +152,7 @@ IF ( BUILD_SERIAL )
 ELSE()
     SET( CTEST_BUILD_COMMAND "${CMAKE_MAKE_PROGRAM} -i -j ${N_PROCS} install" )
 ENDIF()
+SET( CTEST_CUSTOM_WARNING_EXCEPTION "has no symbols" )
 
 
 # Set timeouts: 30 minutes for debug, 15 for opt, and 60 minutes for valgrind/weekly
@@ -158,8 +168,8 @@ ENDIF()
 
 
 # Set valgrind options
-#SET (VALGRIND_COMMAND_OPTIONS "--tool=memcheck --leak-check=yes --track-fds=yes --num-callers=50 --show-reachable=yes --track-origins=yes --malloc-fill=0xff --free-fill=0xfe --suppressions=${LBPM_SOURCE_DIR}/ValgrindSuppresionFile" )
-SET( VALGRIND_COMMAND_OPTIONS  "--tool=memcheck --leak-check=yes --track-fds=yes --num-callers=50 --show-reachable=yes --suppressions=${LBPM_SOURCE_DIR}/ValgrindSuppresionFile" )
+#SET (VALGRIND_COMMAND_OPTIONS "--tool=memcheck --leak-check=yes --track-fds=yes --num-callers=50 --show-reachable=yes --trace-children=yes --track-origins=yes --malloc-fill=0xff --free-fill=0xfe --suppressions=${LBPM_SOURCE_DIR}/ValgrindSuppresionFile" )
+SET( VALGRIND_COMMAND_OPTIONS  "--tool=memcheck --leak-check=yes --num-callers=50 --show-reachable=yes --trace-children=yes --suppressions=${LBPM_SOURCE_DIR}/ValgrindSuppresionFile" )
 IF ( USE_VALGRIND )
     SET( MEMORYCHECK_COMMAND ${VALGRIND_COMMAND} )
     SET( MEMORYCHECKCOMMAND ${VALGRIND_COMMAND} )
@@ -179,16 +189,21 @@ FILE(WRITE "${CTEST_BINARY_DIRECTORY}/CMakeCache.txt" "CTEST_TEST_CTEST:BOOL=1")
 # Set the configure options
 SET( CTEST_OPTIONS )
 SET( CTEST_OPTIONS "-DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE}" )
-SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DCMAKE_C_COMPILER:PATH=${CC};-DCMAKE_C_FLAGS='${C_FLAGS}';" )
-SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DCMAKE_CXX_COMPILER:PATH=${CXX};-DCMAKE_CXX_FLAGS='${CXX_FLAGS}'" )
-SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DMPI_COMPILER:BOOL=true;-DMPIEXEC=${MPIEXEC};-DUSE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=true")
+SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DCMAKE_C_COMPILER:PATH=${CC};-DCMAKE_CXX_COMPILER:PATH=${CXX}" )
+SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DCMAKE_C_FLAGS='${CFLAGS}';-DCMAKE_CXX_FLAGS='${CXXFLAGS}'" )
+SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DLDFLAGS:STRING='${FLAGS}';-DLDLIBS:STRING='${LDLIBS}'" )
+SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DENABLE_GCOV:BOOL=${ENABLE_GCOV}" )
+SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DMPI_COMPILER:BOOL=true;-DMPIEXEC=${MPIEXEC}")
+IF ( NOT USE_VALGRIND )
+    SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DUSE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=true")
+ENDIF()
 IF ( USE_TIMER )
     SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DUSE_TIMER:BOOL=true;-DTIMER_DIRECTORY='${TIMER_DIRECTORY}'" )
 ELSE()
     SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DUSE_TIMER:BOOL=false" )
 ENDIF()
 IF ( USE_CUDA )
-    SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DUSE_CUDA:BOOL=true;-DCUDA_NVCC_FLAGS='${CUDA_FLAGS}';-DCUDA_HOST_COMPILER=${CUDA_HOST_COMPILER}" )
+    SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DUSE_CUDA:BOOL=true;-DCUDA_NVCC_FLAGS='${CUDA_FLAGS}';-DCUDA_HOST_COMPILER=${CUDA_HOST_COMPILER};-DLIB_TYPE=${LIB_TYPE}" )
 ELSE()
     SET( CTEST_OPTIONS "${CTEST_OPTIONS};-DUSE_CUDA:BOOL=false" )
 ENDIF()
@@ -205,15 +220,20 @@ CTEST_CONFIGURE(
     OPTIONS "${CTEST_OPTIONS}"
 )
 CTEST_BUILD()
-IF ( USE_VALGRIND_MATLAB )
-    CTEST_TEST( INCLUDE MATLAB--test_hello_world  PARALLEL_LEVEL ${N_PROCS} )
+IF ( SKIP_TESTS )
+    # Do not run tests
+    SET( CTEST_COVERAGE_COMMAND )
+ELSEIF ( USE_VALGRIND_MATLAB )
+    CTEST_TEST( INCLUDE MATLAB  PARALLEL_LEVEL ${N_PROCS} )
 ELSEIF ( USE_VALGRIND )
-    CTEST_MEMCHECK( EXCLUDE procs   PARALLEL_LEVEL ${N_PROCS} )
+    # CTEST_MEMCHECK( EXCLUDE "(WEEKLY|procs|example--)"  PARALLEL_LEVEL ${N_PROCS} )
+    CTEST_MEMCHECK( EXCLUDE "(WEEKLY|example--)"  PARALLEL_LEVEL ${N_PROCS} )
+ELSEIF ( RUN_WEEKLY )
+    CTEST_TEST( INCLUDE "(WEEKLY|example--)"  PARALLEL_LEVEL ${N_PROCS} )
 ELSE()
-    # CTEST_TEST( EXCLUDE WEEKLY  PARALLEL_LEVEL ${N_PROCS} )
-    CTEST_TEST( PARALLEL_LEVEL ${N_PROCS} )
+    CTEST_TEST( EXCLUDE "(WEEKLY|example--)"  PARALLEL_LEVEL ${N_PROCS} )
 ENDIF()
-IF( CTEST_COVERAGE_COMMAND )
+IF( ENABLE_GCOV )
     CTEST_COVERAGE()
 ENDIF()
 
diff --git a/cmake/libraries.cmake b/cmake/libraries.cmake
index 11ced168..93927533 100644
--- a/cmake/libraries.cmake
+++ b/cmake/libraries.cmake
@@ -1,32 +1,61 @@
 INCLUDE( "${CMAKE_CURRENT_SOURCE_DIR}/cmake/Find_TIMER.cmake" )
+INCLUDE(CheckCCompilerFlag)
+INCLUDE(CheckCXXCompilerFlag)
 
 
-MACRO( CONFIGURE_LINE_COVERAGE )
-    SET ( COVERAGE_LIBS )
+FUNCTION( CONFIGURE_LINE_COVERAGE )
+    SET( COVERAGE_FLAGS )
+    SET( COVERAGE_LIBS )
     IF ( ENABLE_GCOV )
-        ADD_DEFINITIONS ( -fprofile-arcs -ftest-coverage )
-        SET ( COVERAGE_LIBS -lgcov -fprofile-arcs )
-    ENDIF ()
-ENDMACRO()
+        SET( COVERAGE_FLAGS -DUSE_GCOV )
+        SET( CMAKE_REQUIRED_FLAGS ${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage )
+        CHECK_CXX_SOURCE_COMPILES( "int main() { return 0;}" profile-arcs )
+        IF ( profile-arcs )
+            SET( COVERAGE_FLAGS "${COVERAGE_FLAGS} -fprofile-arcs -ftest-coverage" )
+            SET( COVERAGE_LIBS ${COVERAGE_LIBS} -fprofile-arcs )
+        ENDIF()
+        SET( CMAKE_REQUIRED_FLAGS "${CMAKE_CXX_FLAGS} ${COVERAGE_FLAGS} -lgcov" )
+        CHECK_CXX_SOURCE_COMPILES( "int main() { return 0;}" lgcov )
+        IF ( lgcov )
+            SET( COVERAGE_LIBS -lgcov ${COVERAGE_LIBS} )
+        ENDIF()
+        MESSAGE("Enabling coverage:")
+        MESSAGE("   COVERAGE_FLAGS = ${COVERAGE_FLAGS}")
+        MESSAGE("   COVERAGE_LIBS = ${COVERAGE_LIBS}")
+        ADD_DEFINITIONS( ${COVERAGE_FLAGS} )
+        SET( COVERAGE_FLAGS ${COVERAGE_FLAGS} PARENT_SCOPE )
+        SET( COVERAGE_LIBS  ${COVERAGE_LIBS}  PARENT_SCOPE )
+    ENDIF()
+ENDFUNCTION()
 
 
 # Macro to configure CUDA
 MACRO( CONFIGURE_CUDA )
     CHECK_ENABLE_FLAG( USE_CUDA 0 )
     IF( USE_CUDA )
-        SET( CUDA_FLAGS ${CUDA_NVCC_FLAGS} )
-        SET( CUDA_FIND_QUIETLY 1 )
-        INCLUDE ( FindCUDA )
+        # Include FindCUDA
+        INCLUDE( FindCUDA )
         IF ( NOT CUDA_FOUND )
             MESSAGE ( FATAL_ERROR "CUDA not found" )
         ENDIF()
-        SET(CUDA_NVCC_FLAGS ${CUDA_FLAGS})
+        # Initialize the cuda flags
+        IF ( CUDA_FLAGS )
+            SET( CUDA_NVCC_FLAGS "${CUDA_FLAGS}" )
+        ENDIF()
         IF(NOT CUDA_NVCC_FLAGS)
             # Set minimum requirements
             SET( CUDA_NVCC_FLAGS "-arch=sm_20" )
-        ELSE()
-            STRING( REPLACE " " ";" CUDA_NVCC_FLAGS ${CUDA_NVCC_FLAGS} )
         ENDIF()
+        IF( NOT CMAKE_BUILD_TYPE )
+            MESSAGE(FATAL_ERROR "CMAKE_BUILD_TYPE is not set")
+        ELSEIF( ${CMAKE_BUILD_TYPE} STREQUAL "Debug" )
+            SET( CUDA_NVCC_FLAGS "${CUDA_NVCC_FLAGS} -g -O0" )
+        ELSEIF( ${CMAKE_BUILD_TYPE} STREQUAL "Release" )
+            SET( CUDA_NVCC_FLAGS "${CUDA_NVCC_FLAGS} -O3" )
+        ENDIF()
+        SET( CUDA_PROPAGATE_HOST_FLAGS OFF )
+        SET( CUDA_FIND_QUIETLY )
+        STRING( REPLACE " " ";" CUDA_NVCC_FLAGS ${CUDA_NVCC_FLAGS} )
     ENDIF()
     IF ( NOT USE_CUDA )
         MESSAGE( "Not using CUDA" ) 
@@ -47,7 +76,7 @@ ENDMACRO()
 MACRO( CONFIGURE_MIC )
     CHECK_ENABLE_FLAG( USE_MIC 0 )
     ADD_DEFINITIONS ( "-D USE_MIC" ) 
-ENDMACRO() 
+ENDMACRO()
 
 
 # Macro to find and configure the MPI libraries
@@ -72,13 +101,13 @@ MACRO( CONFIGURE_MPI )
                 # User specified the name of the MPI executable
                 SET ( MPIEXEC ${MPI_DIRECTORY}/bin/${MPIEXEC_CMD} )
                 IF ( NOT EXISTS ${MPIEXEC} )
-                    MESSAGE ( FATAL_ERROR "${MPIEXEC_CMD} not found in ${MPI_DIRECTORY}/bin" )
+                    MESSAGE( FATAL_ERROR "${MPIEXEC_CMD} not found in ${MPI_DIRECTORY}/bin" )
                 ENDIF ()
             ELSE ()
                 # Search for the MPI executable in the current directory
                 FIND_PROGRAM ( MPIEXEC  NAMES mpiexec mpirun lamexec  PATHS ${MPI_DIRECTORY}/bin  NO_DEFAULT_PATH )
                 IF ( NOT MPIEXEC )
-                    MESSAGE ( FATAL_ERROR "Could not locate mpi executable" )
+                    MESSAGE( FATAL_ERROR "Could not locate mpi executable" )
                 ENDIF()
             ENDIF ()
             # Set MPI flags
@@ -94,26 +123,33 @@ MACRO( CONFIGURE_MPI )
             # Perform the default search for MPI
             INCLUDE ( FindMPI )
             IF ( NOT MPI_FOUND )
-                MESSAGE ( FATAL_ERROR "Did not find MPI" )
+                MESSAGE( "  MPI_INCLUDE = ${MPI_INCLUDE}" )
+                MESSAGE( "  MPI_LINK_FLAGS = ${MPI_LINK_FLAGS}" )
+                MESSAGE( "  MPI_LIBRARIES = ${MPI_LIBRARIES}" )
+                MESSAGE( FATAL_ERROR "Did not find MPI" )
             ENDIF ()
-            INCLUDE_DIRECTORIES ( ${MPI_INCLUDE_PATH} )
-            SET ( MPI_INCLUDE ${MPI_INCLUDE_PATH} )
+            INCLUDE_DIRECTORIES( "${MPI_INCLUDE_PATH}" )
+            SET( MPI_INCLUDE "${MPI_INCLUDE_PATH}" )
+        ENDIF()
+        # Check if we need to use MPI for serial tests
+        CHECK_ENABLE_FLAG( USE_MPI_FOR_SERIAL_TESTS 0 )
+        # Set defaults if they have not been set
+        IF ( NOT MPIEXEC )
+            SET( MPIEXEC mpirun )
         ENDIF()
-        # Set defaults
         IF ( NOT MPIEXEC_NUMPROC_FLAG )
             SET(MPIEXEC_NUMPROC_FLAG "-n")
         ENDIF()
         # Check if we need to use MPI for serial tests
-        CHECK_ENABLE_FLAG( USE_MPI_FOR_SERIAL_TESTS 0 )
         SET( MPI_CXXFLAGS -DUSE_MPI -I${MPI_INCLUDE} )
         # Set the definitions
-        ADD_DEFINITIONS ( "-D USE_MPI" )  
-        MESSAGE ( "Using MPI" )
-        MESSAGE ( "  MPIEXEC = ${MPIEXEC}" )
-        MESSAGE ( "  MPIEXEC_NUMPROC_FLAG = ${MPIEXEC_NUMPROC_FLAG}" )
-        MESSAGE ( "  MPI_INCLUDE = ${MPI_INCLUDE}" )
-        MESSAGE ( "  MPI_LINK_FLAGS = ${MPI_LINK_FLAGS}" )
-        MESSAGE ( "  MPI_LIBRARIES = ${MPI_LIBRARIES}" )
+        ADD_DEFINITIONS( "-D USE_MPI" )  
+        MESSAGE( "Using MPI" )
+        MESSAGE( "  MPIEXEC = ${MPIEXEC}" )
+        MESSAGE( "  MPIEXEC_NUMPROC_FLAG = ${MPIEXEC_NUMPROC_FLAG}" )
+        MESSAGE( "  MPI_INCLUDE = ${MPI_INCLUDE}" )
+        MESSAGE( "  MPI_LINK_FLAGS = ${MPI_LINK_FLAGS}" )
+        MESSAGE( "  MPI_LIBRARIES = ${MPI_LIBRARIES}" )
     ELSE()
         SET( USE_MPI_FOR_SERIAL_TESTS 0 )
         SET( MPIEXEC "" )
@@ -121,9 +157,9 @@ MACRO( CONFIGURE_MPI )
         SET( MPI_INCLUDE "" )
         SET( MPI_LINK_FLAGS "" )
         SET( MPI_LIBRARIES "" )
-        MESSAGE ( "Not using MPI, all parallel tests will be disabled" )
+        MESSAGE( "Not using MPI, all parallel tests will be disabled" )
     ENDIF()
-ENDMACRO ()
+ENDMACRO()
 
 
 # Macro to configure system-specific libraries and flags
@@ -141,37 +177,46 @@ MACRO( CONFIGURE_SYSTEM )
     # Remove extra library links
     # Get the compiler
     SET_COMPILER ()
-    # Add the static flag if necessary
-    CHECK_ENABLE_FLAG( USE_EXT_STATIC 0 )
-    IF ( USE_EXT_STATIC )
+    CHECK_ENABLE_FLAG( USE_STATIC 0 )
+    IF ( USE_STATIC )
         SET(CMAKE_SHARED_LIBRARY_LINK_C_FLAGS "-static")    # Add static flag
         SET(CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS "-static")  # Add static flag
     ENDIF()
     # Add system dependent flags
-    IF ( USING_MICROSOFT )
+    MESSAGE("System is: ${CMAKE_SYSTEM_NAME}")
+    IF ( ${CMAKE_SYSTEM_NAME} STREQUAL "Windows" )
         # Windows specific system libraries
         SET( SYSTEM_PATHS "C:/Program Files (x86)/Microsoft SDKs/Windows/v7.0A/Lib/x64" 
-                          "C:/Program Files (x86)/Microsoft Visual Studio 8/VC/PlatformSDK/Lib/AMD64" )
-        FIND_LIBRARY ( PSAPI_LIB    NAMES Psapi    PATHS ${SYSTEM_PATHS}  NO_DEFAULT_PATH )
-        FIND_LIBRARY ( DBGHELP_LIB  NAMES DbgHelp  PATHS ${SYSTEM_PATHS}  NO_DEFAULT_PATH )
-        SET( SYSTEM_LIBS ${PSAPI_LIB} ${DBGHELP_LIB} )
+                          "C:/Program Files (x86)/Microsoft Visual Studio 8/VC/PlatformSDK/Lib/AMD64" 
+                          "C:/Program Files (x86)/Microsoft Visual Studio 12.0/Common7/Packages/Debugger/X64" )
+        FIND_LIBRARY( PSAPI_LIB    NAMES Psapi    PATHS ${SYSTEM_PATHS}  NO_DEFAULT_PATH )
+        FIND_LIBRARY( DBGHELP_LIB  NAMES DbgHelp  PATHS ${SYSTEM_PATHS}  NO_DEFAULT_PATH )
+        FIND_LIBRARY( DBGHELP_LIB  NAMES DbgHelp )
+        IF ( PSAPI_LIB ) 
+            ADD_DEFINITIONS( -D PSAPI )
+            SET( SYSTEM_LIBS ${PSAPI_LIB} )
+        ENDIF()
+        IF ( DBGHELP_LIB ) 
+            ADD_DEFINITIONS( -D DBGHELP )
+            SET( SYSTEM_LIBS ${DBGHELP_LIB} )
+        ELSE()
+            MESSAGE( WARNING "Did not find DbgHelp, stack trace will not be availible" )
+        ENDIF()
         MESSAGE("System libs: ${SYSTEM_LIBS}")
     ELSEIF( ${CMAKE_SYSTEM_NAME} STREQUAL "Linux" )
         # Linux specific system libraries
-        CHECK_C_COMPILER_FLAG("-rdynamic" RESULT)
-        IF(RESULT)
-            SET( SYSTEM_LIBS "-lpthread -lz -ldl -rdynamic" )
-        ELSE()
-            SET( SYSTEM_LIBS "-lpthread -lz -ldl" )
+        SET( SYSTEM_LIBS -lz -lpthread -ldl )
+        IF ( NOT USE_STATIC )
+            SET( SYSTEM_LDFLAGS ${SYSTEM_LDFLAGS} -rdynamic )   # Needed for backtrace to print function names
         ENDIF()
         IF ( USING_GCC )
-            SET( SYSTEM_LIBS ${SYSTEM_LIBS} "-lgfortran" )
+            SET( SYSTEM_LIBS ${SYSTEM_LIBS} -lgfortran )
             SET(CMAKE_C_FLAGS   " ${CMAKE_C_FLAGS} -fPIC" )
             SET(CMAKE_CXX_FLAGS " ${CMAKE_CXX_FLAGS} -fPIC" )
         ENDIF()
     ELSEIF( ${CMAKE_SYSTEM_NAME} STREQUAL "Darwin" )
         # Max specific system libraries
-        SET( SYSTEM_LIBS "-lz -ldl" )
+        SET( SYSTEM_LIBS -lz -lpthread -ldl )
     ELSEIF( ${CMAKE_SYSTEM_NAME} STREQUAL "Generic" )
         # Generic system libraries
     ELSE()
@@ -188,4 +233,8 @@ MACRO ( CONFIGURE_LBPM )
     IF ( NOT TEST_MAX_PROCS )
         SET( TEST_MAX_PROCS 32 )
     ENDIF()
+    # Add the correct paths to rpath in case we build shared libraries
+    SET( CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE )
+    SET( CMAKE_BUILD_WITH_INSTALL_RPATH TRUE )
+    SET(CMAKE_INSTALL_RPATH ${CMAKE_INSTALL_RPATH} "${TIMER_DIRECTORY}" "${LBPM_INSTALL_DIR}/lib" )
 ENDMACRO ()
diff --git a/cmake/macros.cmake b/cmake/macros.cmake
index b9528c5b..94bb1350 100644
--- a/cmake/macros.cmake
+++ b/cmake/macros.cmake
@@ -10,7 +10,9 @@ ENDIF()
 IF ( NOT ${PROJ}_INSTALL_DIR )
     MESSAGE(FATAL_ERROR "${PROJ}_INSTALL_DIR must be set before including macros.cmake")
 ENDIF()
-#MESSAGE("Installing project ${PROJ} in ${${PROJ}_INSTALL_DIR}")
+IF ( NOT ${PROJ}_BUILD_DIR )
+    MESSAGE(FATAL_ERROR "${PROJ}_BUILD_DIR must be set before including macros.cmake")
+ENDIF()
 
 
 # Macro to print all variables
@@ -162,37 +164,39 @@ ENDMACRO()
 # Macro to identify the compiler
 MACRO( SET_COMPILER )
     # SET the C/C++ compiler
-    IF( CMAKE_COMPILE_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX )
-        SET( USING_GCC TRUE )
-        ADD_DEFINITIONS( -D USING_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?" )
+    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 )
+            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_MICROSOFT TRUE )
+            ADD_DEFINITIONS( -D USING_MICROSOFT )
+            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 )
+            MESSAGE("Using icc")
+        ELSEIF( ${CMAKE_C_COMPILER_ID} MATCHES "PGI")
+            SET(USING_PGCC TRUE)
+            ADD_DEFINITIONS( -D USING_ICCPGCC )
+            MESSAGE("Using pgCC")
+        ELSEIF( (${CMAKE_C_COMPILER_ID} MATCHES "CRAY") OR (${CMAKE_C_COMPILER_ID} MATCHES "Cray") )
+            SET(USING_CRAY TRUE)
+            ADD_DEFINITIONS( -D USING_CRAY )
+            MESSAGE("Using Cray")
+        ELSEIF( (${CMAKE_C_COMPILER_ID} MATCHES "CLANG") OR (${CMAKE_C_COMPILER_ID} MATCHES "Clang") )
+            SET(USING_CLANG TRUE)
+            ADD_DEFINITIONS( -D USING_CLANG )
+            MESSAGE("Using Clang")
+        ELSE()
+            SET(USING_DEFAULT TRUE)
+            MESSAGE("${CMAKE_C_COMPILER_ID}")
+            MESSAGE("Unknown C/C++ compiler, default flags will be used")
         ENDIF()
-        SET( USING_MICROSOFT TRUE )
-        ADD_DEFINITIONS( -D USING_MICROSOFT )
-        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 )
-        MESSAGE("Using icc")
-    ELSEIF( ${CMAKE_C_COMPILER_ID} MATCHES "PGI")
-        SET(USING_PGCC TRUE)
-        ADD_DEFINITIONS( -D USING_ICCPGCC )
-        MESSAGE("Using pgCC")
-    ELSEIF( (${CMAKE_C_COMPILER_ID} MATCHES "CRAY") OR (${CMAKE_C_COMPILER_ID} MATCHES "Cray") )
-        SET(USING_CRAY TRUE)
-        ADD_DEFINITIONS( -D USING_CRAY )
-        MESSAGE("Using Cray")
-    ELSEIF( (${CMAKE_C_COMPILER_ID} MATCHES "CLANG") OR (${CMAKE_C_COMPILER_ID} MATCHES "Clang") )
-        SET(USING_CLANG TRUE)
-        ADD_DEFINITIONS( -D USING_CLANG )
-        MESSAGE("Using Clang")
-    ELSE()
-        SET(USING_DEFAULT TRUE)
-        MESSAGE("${CMAKE_C_COMPILER_ID}")
-        MESSAGE("Unknown C/C++ compiler, default flags will be used")
     ENDIF()
 ENDMACRO()
 
@@ -202,7 +206,7 @@ MACRO ( SET_COMPILER_FLAGS )
   IF ( USING_GCC )
     # Add gcc specific compiler options
     #    -Wno-reorder:  warning: "" will be initialized after "" when initialized here
-    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+    SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++98")
     SET(CMAKE_C_FLAGS     "${CMAKE_C_FLAGS} -Wall -Wno-char-subscripts -Wno-comment -Wno-unused-variable -Wno-unused-but-set-variable") 
     SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wno-char-subscripts -Wno-comment -Wno-unused-variable -Wno-unused-but-set-variable")
   ELSEIF ( USING_MICROSOFT )
@@ -443,21 +447,40 @@ MACRO( INSTALL_EXAMPLE EXAMPLE )
 ENDMACRO()
 
 
-# Copy an example folder
+# Create an example test
+CONFIGURE_FILE( "${${PROJ}_SOURCE_DIR}/cmake/CompareOutput.cmake" "${${PROJ}_BUILD_DIR}/CompareOutput.cmake" COPYONLY )
 MACRO( TEST_EXAMPLE EXAMPLE EXEFILE PROCS ${ARGN} )
+    SET( EXAMPLE_DIR "${CMAKE_CURRENT_BINARY_DIR}/${EXAMPLE}" )
+    # Copy the example directory
     ADD_CUSTOM_TARGET(
         ${EXAMPLE} ALL
-        ${CMAKE_COMMAND} -E copy_directory "${CMAKE_CURRENT_SOURCE_DIR}/${EXAMPLE}" "${CMAKE_CURRENT_BINARY_DIR}/${EXAMPLE}"
+        ${CMAKE_COMMAND} -E copy_directory "${CMAKE_CURRENT_SOURCE_DIR}/${EXAMPLE}" "${EXAMPLE_DIR}"
         DEPENDS ${EXEFILE}
     )
+    # Create a wrapper script to run the test and copy the output to ${EXAMPLE}.out
+    SET( FILENAME "${EXAMPLE_DIR}/run-${EXAMPLE}" )
+    FILE(WRITE  "${FILENAME}" "# This is a automatically generated file to run example--${EXAMPLE}\n" )
+    FILE(APPEND "${FILENAME}" "${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${PROCS} \"${LBPM_INSTALL_DIR}/bin/${EXEFILE}\" ${ARGN} 2>&1 | tee ${EXAMPLE}.out\n\n" )
+    # Create the test to run the example
     SET( TESTNAME example--${EXAMPLE} )
+    EXECUTE_PROCESS(COMMAND chmod 755 "${FILENAME}")
     ADD_TEST( 
         NAME ${TESTNAME}
-        WORKING_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/${EXAMPLE}" 
-        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${PROCS} "${LBPM_INSTALL_DIR}/bin/${EXEFILE}" ${ARGN} )
+        WORKING_DIRECTORY "${EXAMPLE_DIR}"
+        COMMAND "${FILENAME}" 
+    )
     SET_TESTS_PROPERTIES( ${TESTNAME} PROPERTIES FAIL_REGULAR_EXPRESSION "${TEST_FAIL_REGULAR_EXPRESSION}" PROCESSORS ${PROCS} )
     SET_TESTS_PROPERTIES( ${TESTNAME} PROPERTIES RESOURCE_LOCK ${EXEFILE} )
-
+    # Create a test that checks the output against the data in EXAMPLE/OutputAns.txt
+    IF ( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/${EXAMPLE}/ExampleOutput.txt" )
+        ADD_TEST( 
+            NAME ${TESTNAME}-output
+            WORKING_DIRECTORY "${EXAMPLE_DIR}"
+            COMMAND ${CMAKE_COMMAND} -DTEST=${EXAMPLE}.out -DGOLD=ExampleOutput.txt -P "${${PROJ}_BUILD_DIR}/CompareOutput.cmake"
+        )
+        SET_TESTS_PROPERTIES( ${TESTNAME} PROPERTIES FAIL_REGULAR_EXPRESSION "${TEST_FAIL_REGULAR_EXPRESSION}" PROCESSORS 1 )
+        SET_TESTS_PROPERTIES( ${TESTNAME} PROPERTIES DEPENDS ${TESTNAME} )
+    ENDIF()
 ENDMACRO()
 
 
@@ -539,17 +562,19 @@ MACRO( ADD_DISTCLEAN ${ARGN} )
         Testing
         include
         doc
+        docs
+        latex_docs
         lib
-        tests
+        Makefile.config
+        install_manifest.txt
+        test
+        matlab
+        mex
+        tmp
+        #tmp#
         bin
-        liblbpm-wia.a
-        liblbpm-wia.so
-        cpu
-        gpu
-        example
-        common
-        visit
-        IO
+        cmake
+        ${ARGN}
     )
     ADD_CUSTOM_TARGET (distclean @echo cleaning for source distribution)
     IF (UNIX)
@@ -566,6 +591,7 @@ MACRO( ADD_DISTCLEAN ${ARGN} )
             *.vcxproj*
             ipch
             x64
+            Debug
         )
         SET( DISTCLEAN_FILE "${CMAKE_CURRENT_BINARY_DIR}/distclean.bat" )
         FILE( WRITE  "${DISTCLEAN_FILE}" "del /s /q /f " )
diff --git a/common/Array.h b/common/Array.h
index f3406c02..4e00e81f 100644
--- a/common/Array.h
+++ b/common/Array.h
@@ -181,7 +181,7 @@ IntArray::IntArray(int nx, int ny, int nz)
 
 IntArray::~IntArray()
 {
-	delete data;
+	delete[] data;
 }
 
 
@@ -291,7 +291,7 @@ void DoubleArray::New(int nx, int ny,int nz)
 
 DoubleArray::~DoubleArray()
 {
-	delete data;
+	delete [] data;
 }
 
 double DoubleArray::e(int i)
diff --git a/common/Domain.h b/common/Domain.h
index 17818a6f..07925c66 100755
--- a/common/Domain.h
+++ b/common/Domain.h
@@ -200,6 +200,13 @@ private:
 	}
 };
 
+// Inline function to read line without a return argument
+static inline void fgetl( char * str, int num, FILE * stream )
+{
+    char* ptr = fgets( str, num, stream );
+    if ( 0 ) {char *temp = (char *)&ptr; temp++;}
+}
+
 Domain::~Domain(){
 	delete sendData_x;
 	delete sendData_y;
@@ -775,11 +782,11 @@ inline void ReadSpherePacking(int nspheres, double *List_cx, double *List_cy, do
 	INSIST(fid!=NULL,"Error opening pack.out");
 	//.........Trash the header lines..........
 	char * line = new char[100];
-	fgets(line, 100, fid);
-	fgets(line, 100, fid);
-	fgets(line, 100, fid);
-	fgets(line, 100, fid);
-	fgets(line, 100, fid);
+	fgetl(line, 100, fid);
+	fgetl(line, 100, fid);
+	fgetl(line, 100, fid);
+	fgetl(line, 100, fid);
+	fgetl(line, 100, fid);
 	//........read the spheres..................
     // We will read until a blank like or end-of-file is reached
 	int count = 0;
diff --git a/common/PointList.h b/common/PointList.h
index 66d561d9..2ea39345 100644
--- a/common/PointList.h
+++ b/common/PointList.h
@@ -101,9 +101,9 @@ protected:
 public:
     
     virtual ~DTList() {
-      //  --(*refCount);
-      //  if (*refCount==0) {delete [] Data; delete refCount;}
-      //  Data = 0; refCount = 0; length=0;
+      --(*refCount);
+      if (*refCount==0) {delete [] Data; delete refCount;}
+      Data = 0; refCount = 0; length=0;
     }
 	
     DTList<T> &operator=(const DTList<T> &A) {
diff --git a/common/TwoPhase.h b/common/TwoPhase.h
index f7cc9426..20bc7dcc 100644
--- a/common/TwoPhase.h
+++ b/common/TwoPhase.h
@@ -8,7 +8,6 @@ class TwoPhase{
 	//...........................................................................
 	int n_nw_pts,n_ns_pts,n_ws_pts,n_nws_pts,n_local_sol_pts,n_local_nws_pts;
 	int n_nw_tris,n_ns_tris,n_ws_tris,n_nws_seg,n_local_sol_tris;
-	int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};
 	//...........................................................................
 	int nc;
 	int kstart,kfinish;
@@ -220,7 +219,9 @@ void TwoPhase::ColorToSignedDistance(double Beta, double *ColorData, double *Dis
 	double temp=0.5/Beta;
 	for (int n=0; n<Nx*Ny*Nz; n++){
 		double value = ColorData[n];
-		DistData[n] = temp*log((1.0+value)/(1.0-value));
+		if (value > 0.999 ) DistData[n] = 4.0;
+		else if (value < -0.999 ) DistData[n] = -4.0;
+		else 	DistData[n] = temp*log((1.0+value)/(1.0-value));
 	}
 
 //	for (int n=0; n<Nx*Ny*Nz; n++)	DistData[n] = ColorData[n];
@@ -345,6 +346,8 @@ void TwoPhase::UpdateMeshValues(){
 void TwoPhase::ComputeLocal(){
 	int i,j,k,n;
 	double delphi;
+	int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};
+
 	for (int c=0;c<ncubes;c++){
 		// Get cube from the list
 		i = cubeList(0,c);
diff --git a/common/Utilities.cpp b/common/Utilities.cpp
index 491c8283..0075986a 100644
--- a/common/Utilities.cpp
+++ b/common/Utilities.cpp
@@ -32,12 +32,14 @@
     #include <execinfo.h>
     #include <dlfcn.h>
     #include <mach/mach.h>
+    #include <unistd.h>
 #elif defined(__linux) || defined(__unix) || defined(__posix)
     #define USE_LINUX
     #include <sys/time.h>
     #include <execinfo.h>
     #include <dlfcn.h>
     #include <malloc.h>
+    #include <unistd.h>
 #else
     #error Unknown OS
 #endif
diff --git a/common/pmmc.h b/common/pmmc.h
index 434f9233..aee1b976 100644
--- a/common/pmmc.h
+++ b/common/pmmc.h
@@ -5,6 +5,7 @@
 #include <math.h>
 #include "Array.h"
 #include "PointList.h"
+#include "Utilities.h"
 //#include "vecLib/clapack.h"
 
 using namespace std;
@@ -4159,14 +4160,14 @@ inline void pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, Do
 {
 	int p;
 	double s,lwns,norm;
-	double a,b,c,d,e,f,g,h;
-	double x,y,z,zeta;
+	double zeta;
 	double tangent_x,tangent_y,tangent_z;
 	double ns_x, ns_y, ns_z;
 	double nwn_x, nwn_y, nwn_z;
 	double nwns_x, nwns_y, nwns_z;
 	Point P,A,B;
 	lwns = 0.0;
+    NULL_USE(lwns);
 
 	TriLinPoly Px,Py,Pz,SDx,SDy,SDz,Pt;
 	Px.assign(Fx,i,j,k);
@@ -4259,6 +4260,9 @@ inline void pmmc_CommonCurveSpeed(DoubleArray &CubeValues, DoubleArray &dPdt, Do
 			ReturnVector(2) += zeta*nwns_z*s;
 		}
 	}
+    NULL_USE(tangent_x);
+    NULL_USE(tangent_y);
+    NULL_USE(tangent_z);
 }
 inline void pmmc_CurveOrientation(DoubleArray &Orientation, DTMutableList<Point> &Points, int npts, int i, int j, int k){
 
@@ -4301,10 +4305,9 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, DoubleArray &KN,
 		double &KNavg, double &KGavg, DTMutableList<Point> &Points, int npts, int ic, int jc, int kc){
 	
 	int p,i,j,k;
-	double x,y,z,length;
+	double length;
 	double fxx,fyy,fzz,fxy,fxz,fyz,fx,fy,fz;
 	double sxx,syy,szz,sxy,sxz,syz,sx,sy,sz;
-	
 	double Axx,Axy,Axz,Ayx,Ayy,Ayz,Azx,Azy,Azz;
 //	double Tx[8],Ty[8],Tz[8];	// Tangent vector
 //	double Nx[8],Ny[8],Nz[8];	// Principle normal
@@ -4477,7 +4480,17 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s, DoubleArray &KN,
 			KGavg += K*(nwsx*nwnsx + nwsy*nwnsy + nwsz*nwnsz)*length;
 		}
 	}
+	NULL_USE(fxx); NULL_USE(fyy); NULL_USE(fzz);
+    NULL_USE(fxy); NULL_USE(fxz); NULL_USE(fyz);
+    NULL_USE(fx);  NULL_USE(fy);  NULL_USE(fz);
+	NULL_USE(sxx); NULL_USE(syy); NULL_USE(szz);
+    NULL_USE(sxy); NULL_USE(sxz); NULL_USE(syz);
+    NULL_USE(sx);  NULL_USE(sy);  NULL_USE(sz);
+	NULL_USE(Axx); NULL_USE(Ayy); NULL_USE(Azz);
+    NULL_USE(Axy); NULL_USE(Axz); NULL_USE(Ayz);
 }
+
+
 //--------------------------------------------------------------------------------------------------------
 inline void pmmc_InterfaceSpeed(DoubleArray &dPdt, DoubleArray &P_x, DoubleArray &P_y, DoubleArray &P_z,
 									DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles,
diff --git a/example/Bubble/ExampleOutput.txt b/example/Bubble/ExampleOutput.txt
index 69451121..a80f114f 100644
--- a/example/Bubble/ExampleOutput.txt
+++ b/example/Bubble/ExampleOutput.txt
@@ -1,46 +1,5 @@
-********************************************************
-Running Hybrid Implementation of Color LBM	
-********************************************************
-********************************************************
-tau = 1.000000 
-alpha = 0.010000 
-beta = 0.950000 
-das = 0.100000 
-dbs = 0.900000 
-Value of phi at solid surface = 0.800000 
-Distance to phi = 0.0: 1.156434 
-gamma_{wn} = 0.057960 
-Force(x) = 0.000000 
-Force(y) = 0.000000 
-Force(z) = 0.000000 
-Sub-domain size = 80 x 80 x 80
-Parallel domain size = 1 x 1 x 1
-********************************************************
-Number of blocks = 32 
-Threads per block = 128 
-Sweeps per thread = 135 
-Number of nodes per side = 82 
-Total Number of nodes = 551368 
-********************************************************
-Read input media... 
-Setting up communication control structures 
-Preparing the sendlists 
-SendLists are ready on host
-Prepare to copy send/recv Lists to device 
-Devices are ready to communicate. 
-Copying phase ID to device 
-Allocating distributions 
-********************************************************
-No. of timesteps: 1000 
---------------------------------------------------------------------------------------
-radius sw pw pn awn Jwn Gwn [xx, yy, zz, xy, xz, yz] --------------------------------------------------------------------------------------
-8 0.34906 0.33333 0.0013898 0.26946 0.33334 0.33333 0.33333 -5.3309e-06 1.1201e-05 -1.3209e-05 
-10 0.34568 0.33333 0.0022109 0.21364 0.33334 0.33333 0.33333 -4.1967e-06 1.0265e-05 -1.2941e-05 
-12 0.34355 0.33333 0.0031972 0.1779 0.33334 0.33333 0.33333 -5.6086e-06 1.1842e-05 -1.4445e-05 
-15 0.34143 0.33332 0.005043 0.14164 0.33334 0.33333 0.33333 -5.5737e-06 1.0229e-05 -1.2561e-05 
--------------------------------------------------------------------
-********************************************************
-CPU time = inf 
-Lattice update rate (per core)= 0.000000 MLUPS 
-Lattice update rate (total)= 0.000000 MLUPS 
-********************************************************
+8 0.34914 0.33333 0.001388 0.25364 0.33333 0.33333 0.33333 -1.6662e-06 2.1469e-06 -2.3454e-06
+10 0.34583 0.33333 0.0022065 0.20135 0.33333 0.33333 0.33333 -8.5941e-07 1.734e-06 -2.071e-06
+12 0.34377 0.33332 0.0031884 0.1675 0.33333 0.33333 0.33333 -7.089e-07 1.4714e-06 -1.7527e-06
+15 0.3417 0.33332 0.0050284 0.13345 0.33333 0.33333 0.33333 -7.8226e-07 1.4551e-06 -1.6459e-06
+
diff --git a/sample_scripts/configure_titan_mbt b/sample_scripts/configure_titan_mbt
index 8498f90d..6588a3c0 100755
--- a/sample_scripts/configure_titan_mbt
+++ b/sample_scripts/configure_titan_mbt
@@ -14,9 +14,8 @@ rm -rf CMake*
 cmake                                    \
     -D CMAKE_C_COMPILER:PATH=cc          \
     -D CMAKE_CXX_COMPILER:PATH=CC        \
-    -D CMAKE_CXX_COMPILER:PATH=CC        \
     -D CMAKE_C_FLAGS="-DCBUB"            \
-    -D CMAKE_CXX_FLAGS="-DCBUB"          \
+    -D CMAKE_CXX_FLAGS="-DCBUB -tp=bulldozer"   \
     -D MPI_COMPILER:BOOL=TRUE            \
     -D MPIEXEC=aprun                     \
     -D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
@@ -28,11 +27,11 @@ cmake                                    \
 
 
 # Build the code
-make install -j 8
+# make install -j 8
 
 
 # Run the fast tests
-ctest -E WEEKLY
+# ctest -E WEEKLY
 
 
 # Run the slow tests
diff --git a/tests/BlobIdentify.cpp b/tests/BlobIdentify.cpp
new file mode 100644
index 00000000..7d4f7587
--- /dev/null
+++ b/tests/BlobIdentify.cpp
@@ -0,0 +1,440 @@
+// Sequential blob analysis 
+// Reads parallel simulation data and performs connectivity analysis
+// and averaging on a blob-by-blob basis
+// James E. McClure 2014
+
+#include <iostream>
+#include <math.h>
+#include "pmmc.h"
+//#include "Domain.h"
+
+using namespace std;
+
+inline void ReadCheckpoint(char *FILENAME, double *cDen, double *cDistEven, double *cDistOdd, int N)
+{
+	int q,n;
+	double value;
+	ifstream File(FILENAME,ios::binary);
+	for (n=0; n<N; n++){
+		// Write the two density values
+		File.read((char*) &value, sizeof(value));
+		cDen[n] = value;
+	//	if (n== 66276)	printf("Density a  = %f \n",value);
+		File.read((char*) &value, sizeof(value));
+		cDen[N+n] = value;
+	//	if (n== 66276)	printf("Density b  = %f \n",value);
+		// Read the even distributions
+		for (q=0; q<10; q++){
+			File.read((char*) &value, sizeof(value));
+			cDistEven[q*N+n] = value;
+	//		if (n== 66276)	printf("dist even %i  = %f \n",q,value);
+		}
+		// Read the odd distributions
+		for (q=0; q<9; q++){
+			File.read((char*) &value, sizeof(value));
+			cDistOdd[q*N+n] = value;
+	//		if (n== 66276)	printf("dist even %i  = %f \n",q,value);
+		}
+	}
+	File.close();
+}
+
+inline void ReadBinaryFile(char *FILENAME, double *Data, int N)
+{
+	int n;
+	double value;
+	ifstream File(FILENAME,ios::binary);
+	for (n=0; n<N; n++){
+		// Write the two density values
+		File.read((char*) &value, sizeof(value));
+		Data[n] = value;
+
+	}
+	File.close();
+}
+
+inline void SetPeriodicBC(DoubleArray &Scalar, int nx, int ny, int nz){
+	
+	int i,j,k,in,jn,kn;
+	for (k=0; k<nz; k++){
+		for (j=0; j<ny; j++){
+			for (i=0; i<nx; i++){
+				in = i; jn=j; kn=k;
+				if (i==0) in = nx-2 ;
+				else if (i==nx-1) in = 0;
+				if (j==0) jn = ny-2;
+				else if (j==ny-1) jn = 0;
+				if (k==0) kn = nz-2;
+				else if (k==nz-1) kn = 0;	
+				Scalar(i,j,k) = Scalar(in,jn,kn);
+			}
+		}
+	}
+}
+inline void ReadFromRank(char *FILENAME, DoubleArray &Phase, int nx, int ny, int nz, int iproc, int
+							jproc, int kproc)
+{
+	int i,j,k,q,n,N;
+	int iglobal,jglobal,kglobal;
+	double value;
+	double denA,denB;
+
+	N = nx*ny*nz;
+	
+	double *Den;
+	
+	Den = new double[2*N];
+
+	ifstream File(FILENAME,ios::binary);
+	for (n=0; n<N; n++){
+		// Write the two density values
+		File.read((char*) &value, sizeof(value));
+		Den[2*n] = value;
+		//	if (n== 66276)	printf("Density a  = %f \n",value);
+		File.read((char*) &value, sizeof(value));
+		Den[2*n+1] = value;
+
+		//	if (n== 66276)	printf("Density b  = %f \n",value);
+		// Read the even distributions
+		for (q=0; q<10; q++){
+			File.read((char*) &value, sizeof(value));
+		}
+		// Read the odd distributions
+		for (q=0; q<9; q++){
+			File.read((char*) &value, sizeof(value));
+		}
+	}
+	File.close();
+	
+	// Compute the phase field
+	for (k=1; k<nz-1; k++){
+		for (j=1; j<ny-1; j++){
+			for (i=1; i<nz-1; i++){
+				//........................................................................
+				n = k*nx*ny+j*nx+i;
+				//........................................................................
+				denA = Den[n];
+				denB = Den[N+n];
+				//........................................................................
+				// save values in global arrays
+				//........................................................................
+				iglobal = iproc*(nx-2)+i;
+				jglobal = jproc*(ny-2)+j;
+				kglobal = kproc*(nz-2)+k;
+				//........................................................................				
+				Phase(iglobal,jglobal,kglobal) = (denA-denB)/(denA+denB);
+				//........................................................................
+			}
+		}
+	}
+	
+	delete Den;
+}
+
+int main(int argc, char **argv)
+{
+	printf("-----------------------------------------------------------\n");
+	printf("Labeling Blobs from Two-Phase Lattice Boltzmann Simulation \n");
+	printf("-----------------------------------------------------------\n");
+
+	//.......................................................................
+	int nprocx,nprocy,nprocz,nprocs;
+	int Nx, Ny, Nz;
+	int nx,ny,nz;
+	int nspheres;
+	double Lx,Ly,Lz;
+	//.......................................................................
+	int i,j,k,n,p,idx;
+	int iproc,jproc,kproc;
+	//.......................................................................
+	// Reading the domain information file
+	//.......................................................................
+	ifstream domain("Domain.in");
+	domain >> nprocx;
+	domain >> nprocy;
+	domain >> nprocz;
+	domain >> nx;
+	domain >> ny;
+	domain >> nz;
+	domain >> nspheres;
+	domain >> Lx;
+	domain >> Ly;
+	domain >> Lz;
+	//.......................................................................
+
+	nx+=2;
+	ny+=2;
+	nz+=2;
+	
+	nprocs = nprocx*nprocy*nprocz;
+	printf("Number of MPI ranks: %i \n", nprocs);
+	Nx = (nx-2)*nprocx+2;
+	Ny = (ny-2)*nprocy+2;
+	Nz = (nz-2)*nprocz+2;
+	printf("Full domain size: %i x %i x %i  \n", Nx,Ny,Nz);
+	
+	int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};  // cube corners
+
+	DoubleArray Phase(Nx,Ny,Nz);
+	DoubleArray SignDist(Nx,Ny,Nz);
+	
+	// Filenames used
+	char LocalRankString[8];
+	char LocalRankFilename[40];
+	char LocalRestartFile[40];
+	char BaseFilename[20];
+	char tmpstr[10];
+
+	int proc,iglobal,kglobal,jglobal;
+
+	double * Temp;
+	Temp = new double[nx*ny*nz];
+	
+	for (k=0; k<Nz; k++){
+		for (j=0; j<Ny; j++){
+			for (i=0; i<Nx; i++){
+				SignDist(i,j,k) = -100.0;
+			}
+		}
+	}
+	
+	// read the files and populate main arrays
+	for ( kproc=0; kproc<nprocz; kproc++){
+		for ( jproc=0; jproc<nprocy; jproc++){
+			for ( iproc=0; iproc<nprocx; iproc++){
+				
+				proc = kproc*nprocx*nprocy + jproc*nprocx + iproc;
+
+				sprintf(LocalRankString,"%05d",proc);
+				sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
+				ReadBinaryFile(LocalRankFilename, Temp, nx*ny*nz);	
+				for (k=1; k<nz-1; k++){
+					for (j=1; j<ny-1; j++){
+						for (i=1; i<nz-1; i++){
+
+							//........................................................................
+							n = k*nx*ny+j*nx+i;
+							//........................................................................
+							iglobal = iproc*(nx-2)+i;
+							jglobal = jproc*(ny-2)+j;
+							kglobal = kproc*(nz-2)+k;
+							//........................................................................
+							SignDist(iglobal,jglobal,kglobal) = Temp[n];
+							//........................................................................
+						}
+					}
+				}
+				
+				sprintf(LocalRankFilename,"%s%s","Phase.",LocalRankString);
+				ReadBinaryFile(LocalRankFilename, Temp, nx*ny*nz);	
+				for (k=1; k<nz-1; k++){
+					for (j=1; j<ny-1; j++){
+						for (i=1; i<nx-1; i++){
+
+							//........................................................................
+							n = k*nx*ny+j*nx+i;
+							//........................................................................
+							iglobal = iproc*(nx-2)+i;
+							jglobal = jproc*(ny-2)+j;
+							kglobal = kproc*(nz-2)+k;
+							//........................................................................
+							Phase(iglobal,jglobal,kglobal) = Temp[n];
+							//........................................................................
+						}
+					}
+				}
+				
+				
+			}
+		}
+	}
+	printf("Read %i ranks of Phase, SignDist \n",nprocs);
+	
+	delete [] Temp;
+	
+	IntArray GlobalBlobID(Nx,Ny,Nz);
+
+	SetPeriodicBC(SignDist, Nx, Ny, Nz);
+	SetPeriodicBC(Phase, Nx, Ny, Nz);
+	
+	//		FILE *PHASE;
+	//PHASE = fopen("Phase.dat","wb");
+	//fwrite(Phase.data,8,Nx*Ny*Nz,PHASE);
+	//fclose(PHASE);
+	
+	// Initialize the local blob ID
+	// Initializing the blob ID
+	for (k=0; k<Nz; k++){
+		for (j=0; j<Ny; j++){
+			for (i=0; i<Nx; i++){
+				if (SignDist(i,j,k) < 0.0){
+					// Solid phase 
+					GlobalBlobID(i,j,k) = -2;
+				}
+				else{
+					GlobalBlobID(i,j,k) = -1;
+				}
+			}
+		}
+	}
+	
+	// Compute the porosity
+	double porosity=0.0;
+	for (k=0; k<Nz; k++){
+		for (j=0; j<Ny; j++){
+			for (i=0; i<Nx; i++){
+				if (SignDist(i,j,k) > 0.0){ 
+					porosity += 1.0;
+				}
+			}
+		}
+	}
+	int N=int(porosity*1.25);
+	porosity /= (Nx*Ny*Nz*1.0);
+	printf("Media porosity is %f \n",porosity);
+
+	/* ****************************************************************
+				IDENTIFY ALL BLOBS: F > vF, S > vS
+	****************************************************************** */
+	// Find blob domains, number of blobs
+	int nblobs = 0;					// number of blobs
+	int ncubes = 0;					// total number of nodes in any blob
+	IntArray blobs(3,N);	// store indices for blobs (cubes)
+	IntArray temp(3,N);	// temporary storage array
+	IntArray b(N);          // number of nodes in each blob
+
+	double vF=0.0;
+	double vS=0.0;
+	double trimdist=1.0;
+	printf("Execute blob identification algorithm... \n");
+	// Loop over z=0 first -> blobs attached to this end considered "connected" for LB simulation
+	i=0;
+	int number=0;
+	/*	for (k=0;k<1;k++){
+		for (j=0;j<Ny;j++){
+			if ( Phase(i,j,k) > vF ){
+				if ( SignDist(i,j,k) > vS ){
+					// node i,j,k is in the porespace
+					number = number+ComputeBlob(blobs,nblobs,ncubes,GlobalBlobID,Phase,SignDist,vF,vS,i,j,k,temp);
+				}
+			}
+		}
+	}
+	// Specify the blob on the z axis
+	if (ncubes > 0){
+		b(nblobs) = number;
+//		BlobList.push_back[number];
+		printf("Number of non-wetting phase blobs is: %i \n",nblobs-1);
+		nblobs++;
+	}
+	*/
+	for (k=0;k<Nz;k++){
+		for (j=0;j<Ny;j++){
+			for (i=0;i<Nx;i++){
+				if ( GlobalBlobID(i,j,k) == -1 ){
+					if ( Phase(i,j,k) > vF ){
+						if ( SignDist(i,j,k) > vS ){
+							// node i,j,k is in the porespace
+							b(nblobs) = ComputeBlob(blobs,nblobs,ncubes,GlobalBlobID,Phase,SignDist,vF,vS,i,j,k,temp);
+							nblobs++;							  
+						}
+					}
+				}
+				// Otherwise, this point has already been assigned - ignore
+
+				// Make sure list blob_nodes is large enough
+				if ( nblobs > b.Length-1){
+					printf("Increasing size of blob list \n");
+					b = IncreaseSize(b,b.Length);
+				}
+			}
+		}
+	}
+	// Go over all cubes again -> add any that do not contain nw phase
+	bool add=1;			// Set to false if any corners contain nw-phase ( F > vF)
+//	int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};  // cube corners
+	int count_in=0,count_out=0;
+	int nodx,nody,nodz;
+	for (k=0;k<Nz-1;k++){
+		for (j=0;j<Ny-1;j++){
+			for (i=0;i<Nx-1;i++){
+				// Loop over cube corners
+				add=1;				// initialize to true - add unless corner occupied by nw-phase
+				for (p=0;p<8;p++){
+					nodx=i+cube[p][0];
+					nody=j+cube[p][1];
+					nodz=k+cube[p][2];
+					if ( GlobalBlobID(nodx,nody,nodz) > -1 ){
+						// corner occupied by nw-phase  -> do not add
+						add = 0;
+					}
+				}
+				if ( add == 1 ){
+					blobs(0,ncubes) = i;
+					blobs(1,ncubes) = j;
+					blobs(2,ncubes) = k;
+					ncubes++;
+					count_in++;
+				}
+				else { count_out++; }
+			}
+		}
+	}
+	b(nblobs) = count_in;
+	nblobs++;
+	printf("Identified %i blobs. Writing per-process output files. \n",nblobs);
+
+	int sizeLoc = nx*ny*nz;
+	int *LocalBlobID;
+	LocalBlobID = new int [sizeLoc];
+
+	printf("File size (4 bytes per entry) %i, \n",sizeLoc);
+	// read the files and populate main arrays
+	for ( kproc=0; kproc<nprocz; kproc++){
+		for ( jproc=0; jproc<nprocy; jproc++){
+			for ( iproc=0; iproc<nprocx; iproc++){
+
+				proc = kproc*nprocx*nprocy + jproc*nprocx + iproc;
+
+				sprintf(LocalRankString,"%05d",proc);
+				sprintf(LocalRankFilename,"%s%s","BlobLabel.",LocalRankString);
+
+				for (k=0; k<nz; k++){
+					for (j=0; j<ny; j++){
+						for (i=0; i<nx; i++){
+							//........................................................................
+							n = k*nx*ny+j*nx+i;
+							//........................................................................
+							iglobal = iproc*(nx-2)+i;
+							jglobal = jproc*(ny-2)+j;
+							kglobal = kproc*(nz-2)+k;
+							// periodic BC
+							if (iglobal < 0 ) iglobal+=Nx;
+							if (jglobal < 0 ) jglobal+=Ny;
+							if (kglobal < 0 ) kglobal+=Nz;
+							if (!(iglobal < Nx) ) iglobal-=Nx;
+							if (!(jglobal < Ny) ) jglobal-=Ny;
+							if (!(kglobal < Nz) ) kglobal-=Nz;
+							//........................................................................
+							LocalBlobID[n] = GlobalBlobID(iglobal,jglobal,kglobal);
+							//........................................................................
+						}
+					}
+				}
+
+				FILE *BLOBLOCAL;
+				BLOBLOCAL = fopen(LocalRankFilename,"wb");
+				fwrite(&LocalBlobID[0],4,sizeLoc,BLOBLOCAL);
+				fclose(BLOBLOCAL);
+			}
+		}
+	}
+	printf("Wrote %i ranks of BlobLabel.xxxxx \n",nprocs);
+
+      	FILE *BLOBS;
+	BLOBS = fopen("Blobs.dat","wb");
+	fwrite(GlobalBlobID.data,4,Nx*Ny*Nz,BLOBS);
+	fclose(BLOBS);
+	
+}
+
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 47e56bd7..f1160df5 100755
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -6,6 +6,7 @@ INSTALL_LBPM_EXE( lbpm_disc_pp )
 INSTALL_LBPM_EXE( TestBubble )
 INSTALL_LBPM_EXE( BasicSimulator )
 INSTALL_LBPM_EXE( BlobAnalysis )
+INSTALL_LBPM_EXE( BlobIdentify )
 
 CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_DIR}/cylindertest COPYONLY )
 
diff --git a/tests/TestBubble.cpp b/tests/TestBubble.cpp
index e0c2b25b..0296de47 100644
--- a/tests/TestBubble.cpp
+++ b/tests/TestBubble.cpp
@@ -1155,7 +1155,9 @@ int main(int argc, char **argv)
 	}
 
 	for (int bubbleCount=0; bubbleCount<4; bubbleCount++){
-
+        char bubbleCountName[20];
+        sprintf(bubbleCountName,"bubbleCount-%i",bubbleCount);
+        PROFILE_START(bubbleCountName);
 		BubbleRadius = BubRad[bubbleCount]; // Radius of the current bubble
 
 		// Initialize the bubble
@@ -1356,6 +1358,7 @@ int main(int argc, char **argv)
 		sendtag = recvtag = 5;
 
 		//************ MAIN ITERATION LOOP ***************************************/
+        PROFILE_START("Time-loop");
 		while (timestep < timestepMax){
 
 
@@ -1760,6 +1763,7 @@ int main(int argc, char **argv)
 				WriteCheckpoint(LocalRestartFile, cDen, cDistEven, cDistOdd, N);
 			}
 		}
+        PROFILE_STOP("Time-loop");
 		// End the bubble loop
 		//...........................................................................
 		// Copy the phase indicator field for the later timestep
@@ -2094,7 +2098,7 @@ int main(int argc, char **argv)
 
 
         #ifdef USE_NEW_WRITER
-            std::shared_ptr<IO::TriList> mesh( new IO::TriList() );
+            shared_ptr<IO::TriList> mesh( new IO::TriList() );
             mesh->A.reserve(8*ncubes);
             mesh->B.reserve(8*ncubes);
             mesh->C.reserve(8*ncubes);
@@ -2163,10 +2167,10 @@ int main(int argc, char **argv)
             meshData[0].meshName = "wn-tris";
             meshData[0].mesh = mesh;
             for (size_t k=0; k<meshData.size(); k++) {
-                std::shared_ptr<IO::Variable> dist( new IO::Variable() );
+                shared_ptr<IO::Variable> dist( new IO::Variable() );
                 dist->name = "distance";
                 dist->dim = 1;
-                dist->type = IO::VariableType::NodeVariable;
+                dist->type = IO::NodeVariable;
                 dist->data.resize(3*mesh->A.size());
                 for (size_t i=0; i<mesh->A.size(); i++) {
                     const Point& a = mesh->A[i];
@@ -2183,7 +2187,9 @@ int main(int argc, char **argv)
 		    fclose(WN_TRIS);
         #endif
 
+        PROFILE_STOP(bubbleCountName);
 	}
+
 	//************************************************************************/
 	DeviceBarrier();
 	MPI_Barrier(MPI_COMM_WORLD);
diff --git a/tests/TestInterfaceSpeed.cpp b/tests/TestInterfaceSpeed.cpp
index ba18f0a6..7b5cde5c 100644
--- a/tests/TestInterfaceSpeed.cpp
+++ b/tests/TestInterfaceSpeed.cpp
@@ -1,8 +1,16 @@
 #include <iostream>
 #include <math.h>
-#include "pmmc.h"
-//#include "PointList.h"
-//#include "Array.h"
+
+#include "TwoPhase.h"
+#include "Extras.h"
+#include "D3Q19.h"
+#include "D3Q7.h"
+#include "Color.h"
+#include "common/MPI_Helpers.h"
+#include "Communication.h"
+#include "IO/Mesh.h"
+#include "IO/Writer.h"
+#include "ProfilerApp.h"
 
 #define RADIUS 15
 #define CAPRAD 20
@@ -13,105 +21,38 @@
 
 int main (int argc, char *argv[])
 {
-	//  printf("Radius = %s \n,"RADIUS);  
-	int SIZE = N*N*N;
+	// Initialize MPI
+	int rank,nprocs;
+	MPI_Init(&argc,&argv);
+	MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+	MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
+
+
+	int npx,npy,npz;
+	int i,j,k,n;
 	int Nx,Ny,Nz;
-	Nx = Ny = Nz = N;
-	int i,j,k,p,q,r;
-	
-//	double *Solid; // cylinder
-//	double *Phase; // region of the cylinder	
-//	Solid = new double [SIZE];
-//	Phase = new double [SIZE];
-	DoubleArray SignDist(Nx,Ny,Nz);
-	DoubleArray Phase(Nx,Ny,Nz);
-	DoubleArray Phase_x(Nx,Ny,Nz);
-	DoubleArray Phase_y(Nx,Ny,Nz);
-	DoubleArray Phase_z(Nx,Ny,Nz);
-	DoubleArray Sx(Nx,Ny,Nz);
-	DoubleArray Sy(Nx,Ny,Nz);
-	DoubleArray Sz(Nx,Ny,Nz);
-	DoubleArray GaussCurvature(Nx,Ny,Nz);
-	DoubleArray MeanCurvature(Nx,Ny,Nz);
-	
-	double fluid_isovalue = 0.0;
-	double solid_isovalue = 0.0;
+	double Lx,Ly,Lz;
+	Nx=Ny=Nz=N;
+	npx=npy=npz=1;
+	Lx=Ly=Lz=1.0;
+	int BC=0;	// periodic boundary condition
 
-	/* ****************************************************************
-	 VARIABLES FOR THE PMMC ALGORITHM
-	 ****************************************************************** */
-	//...........................................................................
-	// Averaging variables
-	//...........................................................................
-	double awn,ans,aws,lwns,nwp_volume;
-	double efawns,Jwn;
-	double KNwns,KGwns;
-	double As;
-	double dEs,dAwn,dAns;			 // Global surface energy (calculated by rank=0)
-	double awn_global,ans_global,aws_global,lwns_global,nwp_volume_global;	
-	double As_global;
-	//	bool add=1;			// Set to false if any corners contain nw-phase ( F > fluid_isovalue)
-	int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};  // cube corners
-	//	int count_in=0,count_out=0;
-	//	int nodx,nody,nodz;
-	// initialize lists for vertices for surfaces, common line
-	DTMutableList<Point> nw_pts(20);
-	DTMutableList<Point> ns_pts(20);
-	DTMutableList<Point> ws_pts(20);
-	DTMutableList<Point> nws_pts(20);
-	// initialize triangle lists for surfaces
-	IntArray nw_tris(3,20);
-	IntArray ns_tris(3,20);
-	IntArray ws_tris(3,20);
-	// initialize list for line segments
-	IntArray nws_seg(2,20);
-	
-	DTMutableList<Point> tmp(20);
-	//	IntArray store;
-	
-	int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
-	int n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
-	
-	double s,s1,s2,s3;		// Triangle sides (lengths)
-	Point A,B,C,P;
-	//	double area;
-	
-	// Initialize arrays for local solid surface
-	DTMutableList<Point> local_sol_pts(20);
-	int n_local_sol_pts = 0;
-	IntArray local_sol_tris(3,18);
-	int n_local_sol_tris;
-	DoubleArray values(20);
-	DTMutableList<Point> local_nws_pts(20);
-	int n_local_nws_pts;
-	
-	DoubleArray CubeValues(2,2,2);
-	DoubleArray ContactAngle(20);
-	DoubleArray KGwns_values(20);
-	DoubleArray KNwns_values(20);
-	DoubleArray wn_curvature(20);
-	DoubleArray InterfaceSpeed(20);
-	DoubleArray NormalVector(60);
-	DoubleArray vawn(6);
-	DoubleArray vawns(3);
-	
-	for (i=0;i<3;i++){
-		vawn(i) = 0.0;
-		vawns(i) = 0.0;
-	}
+	Domain Dm(Nx,Ny,Nz,rank,npx,npy,npz,Lx,Ly,Lz,BC);
+
+	for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
+
+	Dm.CommInit(MPI_COMM_WORLD);
+
+	TwoPhase Averages(Dm);
+	int timestep=0;
+
+	Nx = Dm.Nx;
+	Ny = Dm.Ny;
+	Nz = Dm.Nz;
 
-	int c;
-	//...........................................................................
-	int ncubes = (Nx-2)*(Ny-2)*(Nz-2);	// Exclude the "upper" halo
-	IntArray cubeList(3,ncubes);
-	pmmc_CubeListFromMesh(cubeList, ncubes, Nx, Ny, Nz);
-	//...........................................................................
 	double Cx,Cy,Cz;
 	double dist1,dist2;
-	// Extra copies of phase indicator needed to compute time derivatives on CPU
-	DoubleArray Phase_tminus(Nx,Ny,Nz);
-	DoubleArray Phase_tplus(Nx,Ny,Nz);
-	DoubleArray dPdt(Nx,Ny,Nz);
+
 	Cx = Cy = Cz = N*0.5;
 	for (k=0; k<Nz; k++){
 		for (j=0; j<Ny; j++){
@@ -119,7 +60,7 @@ int main (int argc, char *argv[])
 				dist2 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)+(k-Cz)*(k-Cz)) - CAPRAD;
 				dist2 = fabs(Cz-k)-HEIGHT;
 
-				Phase_tminus(i,j,k) = dist2;
+				Averages.Phase_tminus(i,j,k) = dist2;
 			}
 		} 
 	}
@@ -132,8 +73,9 @@ int main (int argc, char *argv[])
 				dist2 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)+(k-Cz)*(k-Cz)) - CAPRAD;
 				dist2 = fabs(Cz-k)-HEIGHT;
 
-				SignDist(i,j,k) = -dist1;
-				Phase(i,j,k) = dist2;
+				Averages.SDs(i,j,k) = -dist1;
+				Averages.Phase(i,j,k) = dist2;
+				Averages.SDn(i,j,k) = dist2;
 			}
 		}   
 	}
@@ -144,106 +86,57 @@ int main (int argc, char *argv[])
 				dist2 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)+(k-Cz)*(k-Cz)) - CAPRAD;
 				dist2 = fabs(Cz-k)-HEIGHT;
 
-				Phase_tplus(i,j,k) = dist2;
+				Averages.Phase_tplus(i,j,k) = dist2;
 			}
 		}   
 	}
 	
 	//...........................................................................
-	// Calculate the time derivative of the phase indicator field
-	for (int n=0; n<Nx*Ny*Nz; n++)	dPdt(n) = 0.5*(Phase_tplus(n) - Phase_tminus(n));
-	
-	pmmc_MeshGradient(Phase,Phase_x,Phase_y,Phase_z,Nx,Ny,Nz);
-	pmmc_MeshGradient(SignDist,Sx,Sy,Sz,Nx,Ny,Nz);
-	
-	double norm;
-	for (k=0; k<Nz; k++){
-		for (j=0; j<Ny; j++){
-			for (i=0; i<Nx; i++){
-				norm = Phase_x(i,j,k)*Phase_x(i,j,k)+Phase_y(i,j,k)*Phase_y(i,j,k)+Phase_z(i,j,k)*Phase_z(i,j,k);
-			}
-		}
-	}
-	
-	// End of the loop to set the values
-	awn = aws = ans = lwns = 0.0;
-	nwp_volume = 0.0;
-	As = 0.0;
-	
-	for (c=0;c<ncubes;c++){
-		// Get cube from the list
-		i = cubeList(0,c);
-		j = cubeList(1,c);
-		k = cubeList(2,c);
-		
-		for (p=0;p<8;p++){
-			if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 
-				&&  SignDist(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
-				nwp_volume += 0.125;
-			}
-		}
-		
-		// Run PMMC
-		n_local_sol_tris = 0;
-		n_local_sol_pts = 0;
-		n_local_nws_pts = 0;
-		
-		n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
-		n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
 
-		// Construct the interfaces and common curve
-		pmmc_ConstructLocalCube(SignDist, Phase, solid_isovalue, fluid_isovalue,
-				nw_pts, nw_tris, values, ns_pts, ns_tris, ws_pts, ws_tris,
-				local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
-				n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
-				n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
-				i, j, k, Nx, Ny, Nz);
-
-		efawns += pmmc_CubeContactAngle(CubeValues,ContactAngle,Phase_x,Phase_y,Phase_z,Sx,Sy,Sz,local_nws_pts,i,j,k,n_local_nws_pts);
-		
-		Jwn += pmmc_CubeSurfaceInterpValue(CubeValues, MeanCurvature, nw_pts, nw_tris,
-									wn_curvature, i, j, k, n_nw_pts, n_nw_tris);
-		
-		pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris,
-							NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
-		
-		pmmc_CommonCurveSpeed(CubeValues, dPdt, vawns, Phase_x,Phase_y,Phase_z,Sx,Sy,Sz,
-				local_nws_pts,i,j,k,n_local_nws_pts);
-
-		pmmc_CurveCurvature(Phase, SignDist, KNwns_values, KGwns_values, KNwns, KGwns, nws_pts, n_nws_pts, i, j, k);	
-		
-	//	if (n_nw_pts>0) printf("speed %f \n",InterfaceSpeed(0));
-		
-		//*******************************************************************
-		// Compute the Interfacial Areas, Common Line length
-		awn += pmmc_CubeSurfaceArea(nw_pts,nw_tris,n_nw_tris);
-		ans += pmmc_CubeSurfaceArea(ns_pts,ns_tris,n_ns_tris);
-		aws += pmmc_CubeSurfaceArea(ws_pts,ws_tris,n_ws_tris);
-		As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris);
-		lwns +=  pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
-	}
-	KGwns /= lwns;
-	KNwns /= lwns;
-	Jwn /= awn;
-	efawns /= lwns;
-	for (i=0;i<6;i++)	vawn(i) /= awn;
-	for (i=0;i<3;i++)	vawns(i) /= lwns;
+	//....................................................................
+	// The following only need to be done once
+	Averages.SetupCubes(Dm);
+	Averages.UpdateSolid(); 	// unless the solid is deformable!
+	//....................................................................
+	// The following need to be called each time new averages are computed
+	Averages.Initialize();
+	Averages.UpdateMeshValues();
+	Averages.ComputeLocal();
+	Averages.Reduce();
+	Averages.PrintAll(timestep);
+	//....................................................................
 	
 	printf("-------------------------------- \n");
-	printf("NWP volume = %f \n", nwp_volume);
-	printf("Area wn = %f, Analytical = %f \n", awn,2*PI*RADIUS*RADIUS);
-	printf("Area ns = %f, Analytical = %f \n", ans, 2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT);
-	printf("Area ws = %f, Analytical = %f \n", aws, 4*PI*RADIUS*HEIGHT);
-	printf("Area s = %f, Analytical = %f \n", As, 2*PI*RADIUS*(N-2));
-	printf("Length wns = %f, Analytical = %f \n", lwns, 4*PI*RADIUS);
-	printf("Geodesic curvature (wns) = %f, Analytical = %f \n", KGwns, 0.0);
-	printf("Normal curvature (wns) = %f, Analytical = %f \n", KNwns, 1.0/RADIUS);
+	printf("NWP volume = %f \n", Averages.nwp_volume);
+	printf("Area wn = %f, Analytical = %f \n", Averages.awn,2*PI*RADIUS*RADIUS);
+	printf("Area ns = %f, Analytical = %f \n", Averages.ans, 2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT);
+	printf("Area ws = %f, Analytical = %f \n", Averages.aws, 4*PI*RADIUS*HEIGHT);
+	printf("Area s = %f, Analytical = %f \n", Averages.As, 2*PI*RADIUS*(N-2));
+	printf("Length wns = %f, Analytical = %f \n", Averages.lwns, 4*PI*RADIUS);
+	printf("Geodesic curvature (wns) = %f, Analytical = %f \n", Averages.KGwns_global, 0.0);
+	printf("Normal curvature (wns) = %f, Analytical = %f \n", Averages.KNwns_global, 1.0/RADIUS);
 //	printf("Cos(theta_wns) = %f, Analytical = %f \n",efawns/lwns,1.0*RADIUS/CAPRAD);
-	printf("Interface Velocity = %f,%f,%f \n",vawn(0),vawn(1),vawn(2));
-	printf("Common Curve Velocity = %f,%f,%f \n",vawns(0),vawns(1),vawns(2));
+	printf("Interface Velocity = %f,%f,%f \n",Averages.vawn_global(0),Averages.vawn_global(1),Averages.vawn_global(2));
+	printf("Common Curve Velocity = %f,%f,%f \n",Averages.vawns_global(0),Averages.vawns_global(1),Averages.vawns_global(2));
 	printf("-------------------------------- \n");	
 	//.........................................................................	
 	
 	int toReturn = 0;
+
+	if ( fabs(Averages.vawn_global(2)+0.2) > 0.01){
+		printf("TestInterfaceSpeed: Error too high for kinematic velocity of wn interface \n");
+		toReturn=1;
+	}
+	if ( fabs(Averages.vawns_global(2)+0.2) > 0.01){
+		printf("TestInterfaceSpeed: Error too high for kinematic velocity of common curve \n");
+		toReturn=2;
+	}
+
 	return toReturn;
+
+	// ****************************************************
+	MPI_Barrier(MPI_COMM_WORLD);
+	return 0;
+	MPI_Finalize();
+	// ****************************************************
 }
diff --git a/tests/TestWriter.cpp b/tests/TestWriter.cpp
index dd60e5a8..bb2bc6fc 100644
--- a/tests/TestWriter.cpp
+++ b/tests/TestWriter.cpp
@@ -4,8 +4,8 @@
 #include <exception>
 #include <stdexcept>
 #include <fstream>
-#include <memory>
 
+#include "shared_ptr.h"
 #include "common/UnitTest.h"
 #include "common/Utilities.h"
 #include "common/MPI_Helpers.h"
@@ -55,19 +55,19 @@ int main(int argc, char **argv)
     };
 
     // Create the meshes
-    std::shared_ptr<IO::PointList> set1( new IO::PointList(N_points) );
+    shared_ptr<IO::PointList> set1( new IO::PointList(N_points) );
     for (int i=0; i<N_points; i++) {
         set1->points[i].x = x[i];
         set1->points[i].y = y[i];
         set1->points[i].z = z[i];
     }
-    std::shared_ptr<IO::TriMesh> trimesh( new IO::TriMesh(N_tri,set1) );
+    shared_ptr<IO::TriMesh> trimesh( new IO::TriMesh(N_tri,set1) );
     for (int i=0; i<N_tri; i++) {
         trimesh->A[i] = tri[i][0];
         trimesh->B[i] = tri[i][1];
         trimesh->C[i] = tri[i][2];
     }
-    std::shared_ptr<IO::TriList> trilist( new IO::TriList(*trimesh) );
+    shared_ptr<IO::TriList> trilist( new IO::TriList(*trimesh) );
     for (int i=0; i<N_tri; i++) {
         Point A(x[tri[i][0]],y[tri[i][0]],z[tri[i][0]]);
         Point B(x[tri[i][1]],y[tri[i][1]],z[tri[i][1]]);
@@ -80,14 +80,14 @@ int main(int argc, char **argv)
     }
 
     // Create the variables
-    std::shared_ptr<IO::Variable> dist_set1( new IO::Variable() );
-    std::shared_ptr<IO::Variable> dist_list( new IO::Variable() );
+    shared_ptr<IO::Variable> dist_set1( new IO::Variable() );
+    shared_ptr<IO::Variable> dist_list( new IO::Variable() );
     dist_set1->dim = 1;
     dist_list->dim = 1;
     dist_set1->name = "Distance";
     dist_list->name = "Distance";
-    dist_set1->type = IO::VariableType::NodeVariable;
-    dist_list->type = IO::VariableType::NodeVariable;
+    dist_set1->type = IO::NodeVariable;
+    dist_list->type = IO::NodeVariable;
     dist_set1->data.resize( N_points );
     for (int i=0; i<N_points; i++)
         dist_set1->data[i] = distance(set1->points[i]);
@@ -143,16 +143,16 @@ int main(int argc, char **argv)
         // For each domain, load the mesh and check its data
         for (size_t j=0; j<list.size(); j++) {
             for (size_t k=0; k<list[i].domains.size(); k++) {
-                std::shared_ptr<IO::Mesh> mesh = IO::getMesh(".",timesteps[i],list[j],k);
-                if ( mesh==NULL ) {
+                shared_ptr<IO::Mesh> mesh = IO::getMesh(".",timesteps[i],list[j],k);
+                if ( mesh.get()==NULL ) {
                     printf("Failed to load %s\n",meshData[i].meshName.c_str());
                     pass = false;
                     break;
                 }
                 if ( meshData[j].meshName=="pointmesh" ) {
                     // Check the pointmesh
-                    std::shared_ptr<IO::PointList> pmesh = IO::getPointList(mesh);
-                    if ( pmesh==NULL ) {
+                    shared_ptr<IO::PointList> pmesh = IO::getPointList(mesh);
+                    if ( pmesh.get()==NULL ) {
                         pass = false;
                         break;
                     }
@@ -163,9 +163,9 @@ int main(int argc, char **argv)
                 }
                 if ( meshData[j].meshName=="trimesh" || meshData[j].meshName=="trilist" ) {
                     // Check the trimesh/trilist
-                    std::shared_ptr<IO::TriMesh> mesh1 = IO::getTriMesh(mesh);
-                    std::shared_ptr<IO::TriList> mesh2 = IO::getTriList(mesh);
-                    if ( mesh1==NULL || mesh2==NULL ) {
+                    shared_ptr<IO::TriMesh> mesh1 = IO::getTriMesh(mesh);
+                    shared_ptr<IO::TriList> mesh2 = IO::getTriList(mesh);
+                    if ( mesh1.get()==NULL || mesh2.get()==NULL ) {
                         pass = false;
                         break;
                     }
@@ -213,7 +213,7 @@ int main(int argc, char **argv)
             }
             for (size_t v=0; v<list[i].variables.size(); v++) {
                 for (size_t k=0; k<list[i].domains.size(); k++) {
-                    std::shared_ptr<const IO::Variable> variable = 
+                    shared_ptr<const IO::Variable> variable = 
                         IO::getVariable(".",timesteps[i],list[j],k,list[j].variables[v].name);
                     const IO::Variable& var1 = *mesh0->vars[v];
                     const IO::Variable& var2 = *variable;
diff --git a/tests/hello_world.cpp b/tests/hello_world.cpp
index 6fce2a48..b67e7e6c 100644
--- a/tests/hello_world.cpp
+++ b/tests/hello_world.cpp
@@ -1,5 +1,6 @@
 #include <iostream>
 #include "common/MPI_Helpers.h"
+#include "common/Utilities.h"
 
 
 int main (int argc, char **argv)
@@ -15,6 +16,12 @@ int main (int argc, char **argv)
         MPI_Barrier(MPI_COMM_WORLD);
     }
 
+    // Create a memory leak for valgrind to find
+    if ( nprocs==1 ) {
+        double *x = new double[1];
+        ASSERT(x!=NULL);
+    }
+
     // set the error code
     // Note: the error code should be consistent across all processors
     int error = 0;
diff --git a/tests/lb2_Color_wia_mpi.cpp b/tests/lb2_Color_wia_mpi.cpp
index 2562fcc0..47cd342d 100644
--- a/tests/lb2_Color_wia_mpi.cpp
+++ b/tests/lb2_Color_wia_mpi.cpp
@@ -509,7 +509,7 @@ int main(int argc, char **argv)
 #endif
 
 	// Initialize communication structures in averaging domain
-	for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
+	for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = id[i];
 	Dm.CommInit(MPI_COMM_WORLD);
 
 	// Set up MPI communication structures
diff --git a/tests/lbpm_BlobAnalysis.cpp b/tests/lbpm_BlobAnalysis.cpp
index 766b4a8f..772274c9 100644
--- a/tests/lbpm_BlobAnalysis.cpp
+++ b/tests/lbpm_BlobAnalysis.cpp
@@ -144,7 +144,7 @@ int main(int argc, char **argv)
 	//.................................................
 	Domain Dm(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC);
 	TwoPhase Averages(Dm);
-    BlobTwoPhase BlobAverages(nblobs_global);
+	BlobTwoPhase BlobAverages(nblobs_global);
 	//.......................................................................
 	// Filenames used
 	char LocalRankString[8];
@@ -223,7 +223,6 @@ int main(int argc, char **argv)
 	for (k=0;k<Nz;k++){
 		for (j=0;j<Ny;j++){
 			for (i=1;i<Nx;i++){
-                
 			}
 		}
 	}
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index a920182a..2c39098d 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -424,7 +424,7 @@ int main(int argc, char **argv)
 		FlipID(id,Nx*Ny*Nz);
 	}
 	// Initialize communication structures in averaging domain
-	for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = 1;
+	for (i=0; i<Dm.Nx*Dm.Ny*Dm.Nz; i++) Dm.id[i] = id[i];
 	Dm.CommInit(MPI_COMM_WORLD);
 
 	// Set up MPI communication structurese
@@ -1547,7 +1547,7 @@ int main(int argc, char **argv)
 		// Timestep completed!
 		timestep++;
 		//...................................................................
-		if (timestep%1000 == 995){
+		if (timestep%5000 == 995){
 			//...........................................................................
 			// Copy the phase indicator field for the earlier timestep
 			DeviceBarrier();
@@ -1555,7 +1555,7 @@ int main(int argc, char **argv)
 			Averages.ColorToSignedDistance(beta,Averages.Phase.data,Averages.Phase_tplus.data);
 			//...........................................................................
 		}
-		if (timestep%1000 == 0){
+		if (timestep%5000 == 0){
 			//...........................................................................
 			// Copy the data for for the analysis timestep
 			//...........................................................................
@@ -1570,7 +1570,7 @@ int main(int argc, char **argv)
 			CopyToHost(Averages.Vel_z.data,&Velocity[2*N],N*sizeof(double));
 			MPI_Barrier(MPI_COMM_WORLD);
 		}
-		if (timestep%1000 == 5){
+		if (timestep%5000 == 5){
 			//...........................................................................
 			// Copy the phase indicator field for the later timestep
 			DeviceBarrier();
diff --git a/tests/testUtilities.cpp b/tests/testUtilities.cpp
index 09a516d3..99ce09d6 100644
--- a/tests/testUtilities.cpp
+++ b/tests/testUtilities.cpp
@@ -1,18 +1,31 @@
-#include <stdlib.h>
-#include <stdio.h>
 #include <iostream>
+#include <vector>
+#include <stdio.h>
+#include <stdlib.h>
 #include <time.h>
 #include <sys/stat.h>
 #include <math.h>
 #include <stdexcept>
-#include <vector>
-#include <string>
+#include <string.h>
+#include <stdint.h>
 
 #include "common/Utilities.h"
 #include "common/UnitTest.h"
 #include "common/MPI_Helpers.h"
 
 
+// Detect the OS (defines which tests we allow to fail)
+#if defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64) || defined(_MSC_VER)
+    #define USE_WINDOWS
+#elif defined(__APPLE__)
+    #define USE_MAC
+#elif defined(__linux) || defined(__unix) || defined(__posix)
+    #define USE_LINUX
+#else
+    #error Unknown OS
+#endif
+
+
 // Function to return the call stack
 std::vector<std::string> get_call_stack() 
 {
@@ -32,73 +45,88 @@ int main(int argc, char *argv[])
     UnitTest ut;
     Utilities::setAbortBehavior( true, true, true );
 
-    // Test the memory usage
-    double t0 = Utilities::time();
-    size_t n_bytes1 = Utilities::getMemoryUsage();
-    double time1 = Utilities::time() - t0;
-    double *tmp = new double[0x100000];
-    NULL_USE(tmp);
-    t0 = Utilities::time();
-    size_t n_bytes2 = Utilities::getMemoryUsage();
-    double time2 = Utilities::time() - t0;
-    delete [] tmp;
-    t0 = Utilities::time();
-    size_t n_bytes3 = Utilities::getMemoryUsage();
-    double time3 = Utilities::time() - t0;
-    std::cout << "Number of bytes used for a basic test: " << n_bytes1 << ", " << n_bytes2 << ", " << n_bytes3 << std::endl;
-    std::cout << "   Time to query: " << time1*1e6 << " us, " << time2*1e6 << " us, " << time3*1e6 << " us" << std::endl;
-    if ( n_bytes1==0 ) {
-        ut.failure("getMemoryUsage returns 0");
-    } else {
-        ut.passes("getMemoryUsage returns non-zero");
-        if ( n_bytes2>n_bytes1 )
-            ut.passes("getMemoryUsage increases size");
-        else
-            ut.failure("getMemoryUsage increases size");
-        if ( n_bytes1==n_bytes3 )
-            ut.passes("getMemoryUsage decreases size properly");
-        else
-            ut.expected_failure("getMemoryUsage does not decrease size properly");
-    }
+    // Limit the scope of variables
+    { 
+        // Test the memory usage
+        double t0 = Utilities::time();
+        size_t n_bytes1 = Utilities::getMemoryUsage();
+        double time1 = Utilities::time() - t0;
+        uint64_t *tmp = new uint64_t[0x100000];
+        memset(tmp,0xAA,0x100000*sizeof(uint64_t));
+        NULL_USE(tmp);
+        t0 = Utilities::time();
+        size_t n_bytes2 = Utilities::getMemoryUsage();
+        double time2 = Utilities::time() - t0;
+        delete [] tmp;
+        t0 = Utilities::time();
+        size_t n_bytes3 = Utilities::getMemoryUsage();
+        double time3 = Utilities::time() - t0;
+        std::cout << "Number of bytes used for a basic test: " << n_bytes1 << ", " << n_bytes2 << ", " << n_bytes3 << std::endl;
+        std::cout << "   Time to query: " << time1*1e6 << " us, " << time2*1e6 << " us, " << time3*1e6 << " us" << std::endl;
+        if ( n_bytes1==0 ) {
+            ut.failure("getMemoryUsage returns 0");
+        } else {
+            ut.passes("getMemoryUsage returns non-zero");
+            if ( n_bytes2>n_bytes1 ) {
+                ut.passes("getMemoryUsage increases size");
+            } else {
+                #if defined(USE_MAC)
+                    ut.expected_failure("getMemoryUsage does not increase size");
+                #else
+                    ut.failure("getMemoryUsage increases size");
+                #endif
+            }
+            if ( n_bytes1==n_bytes3 ) {
+                ut.passes("getMemoryUsage decreases size properly");
+            } else {
+                #if defined(USE_MAC) || defined(USE_WINDOWS)
+                    ut.expected_failure("getMemoryUsage does not decrease size properly");
+                #else
+                    ut.failure("getMemoryUsage does not decrease size properly");
+                #endif
+            }
+        }
 
-    // Test getting the current call stack
-    std::vector<std::string> call_stack = get_call_stack();
-    if ( rank==0 ) {
-        std::cout << "Call stack:" << std::endl;
-        for (size_t i=0; i<call_stack.size(); i++)
-            std::cout << "   " << call_stack[i] << std::endl;
-    }
-    if ( !call_stack.empty() ) {
-        ut.passes("non empty call stack");
-        if ( call_stack[0].find("get_call_stack()") != std::string::npos )
-            ut.passes("call stack decoded function symbols");
-        else
-            ut.expected_failure("call stack was unable to decode function symbols");
-    } else {
-        ut.failure("non empty call stack");
-    }
+        // Test getting the current call stack
+        std::vector<std::string> call_stack = get_call_stack();
+        if ( rank==0 ) {
+            std::cout << "Call stack:" << std::endl;
+            for (size_t i=0; i<call_stack.size(); i++)
+                std::cout << "   " << call_stack[i] << std::endl;
+        }
+        if ( !call_stack.empty() ) {
+            ut.passes("non empty call stack");
+            if ( call_stack[0].find("get_call_stack()") != std::string::npos )
+                ut.passes("call stack decoded function symbols");
+            else
+                ut.expected_failure("call stack was unable to decode function symbols");
+        } else {
+            ut.failure("non empty call stack");
+        }
+
+        // Test catching an error
+        try {
+            ERROR("Test");
+            ut.failure("Failed to catch RAY_ERROR");
+        } catch (...) {
+            ut.passes("Caught RAY_ERROR");
+        }
+        try {
+            throw std::logic_error("test");
+            ut.failure("Failed to catch exception");
+        } catch (...) {
+            ut.passes("Caught exception");
+        }
+        
+        // Test time/tick
+        double time = Utilities::time();
+        double res = Utilities::tick();
+        if ( time==0 || res==0 )
+            ut.failure("time/tick");
+        else
+            ut.passes("time/tick");
 
-    // Test catching an error
-    try {
-        ERROR("Test");
-        ut.failure("Failed to catch RAY_ERROR");
-    } catch (...) {
-        ut.passes("Caught RAY_ERROR");
     }
-    try {
-        throw std::logic_error("test");
-        ut.failure("Failed to catch exception");
-    } catch (...) {
-        ut.passes("Caught exception");
-    }
-    
-    // Test time/tick
-    double time = Utilities::time();
-    double res = Utilities::tick();
-    if ( time==0 || res==0 )
-        ut.failure("time/tick");
-    else
-        ut.passes("time/tick");
 
     // Finished
     ut.report();