diff --git a/analysis/dcel.cpp b/analysis/dcel.cpp
index 01e7b571..71f69951 100644
--- a/analysis/dcel.cpp
+++ b/analysis/dcel.cpp
@@ -1,7 +1,5 @@
 #include "analysis/dcel.h"
 
-
-
 DECL::DECL(){
 }
 
@@ -369,243 +367,37 @@ double DECL::EdgeAngle(int edge)
 	return angle;
 }
 
-void Isosurface(DoubleArray &A, const double &v)
+void iso_surface(const Array<double>&Field, const double isovalue)
 {
-    NULL_USE( v );
-
-	Point P,Q;
-	Point PlaceHolder;
-	Point C0,C1,C2,C3,C4,C5,C6,C7;
-
-	int TriangleCount;
-	int VertexCount;
-	int CubeIndex;
-
-	Point VertexList[12];
-	Point NewVertexList[12];
-	int LocalRemap[12];
-
-	Point cellvertices[20];
-    std::array<std::array<int,3>,20> Triangles;
-    Triangles.fill( { 0 } );
-
-	// Values from array 'A' at the cube corners
-	double CubeValues[8];
-
-	int Nx = A.size(0);
-	int Ny = A.size(1);
-	int Nz = A.size(2);
-
-	// Points corresponding to cube corners
-	C0.x = 0.0; C0.y = 0.0; C0.z = 0.0;
-	C1.x = 1.0; C1.y = 0.0; C1.z = 0.0;
-	C2.x = 1.0; C2.y = 1.0; C2.z = 0.0;
-	C3.x = 0.0; C3.y = 1.0; C3.z = 0.0;
-	C4.x = 0.0; C4.y = 0.0; C4.z = 1.0;
-	C5.x = 1.0; C5.y = 0.0; C5.z = 1.0;
-	C6.x = 1.0; C6.y = 1.0; C6.z = 1.0;
-	C7.x = 0.0; C7.y = 1.0; C7.z = 1.0;							
-
-	std::vector<std::array<int,6>> HalfEdge;
-	for (int k=1; k<Nz-1; k++){		
-		for (int j=1; j<Ny-1; j++){
-			for (int i=1; i<Nx-1; i++){
-				// Set the corner values for this cube
-				CubeValues[0] = A(i,j,k);
-				CubeValues[1] = A(i+1,j,k);
-				CubeValues[2] = A(i+1,j+1,k);
-				CubeValues[3] = A(i,j+1,k);
-				CubeValues[4] = A(i,j,k+1);
-				CubeValues[5] = A(i+1,j,k+1);
-				CubeValues[6] = A(i+1,j+1,k+1);
-				CubeValues[7] = A(i,j+1,k+1);
-
-				//Determine the index into the edge table which
-				//tells us which vertices are inside of the surface
-				CubeIndex = 0;
-				if (CubeValues[0] < 0.0f) CubeIndex |= 1;
-				if (CubeValues[1] < 0.0f) CubeIndex |= 2;
-				if (CubeValues[2] < 0.0f) CubeIndex |= 4;
-				if (CubeValues[3] < 0.0f) CubeIndex |= 8;
-				if (CubeValues[4] < 0.0f) CubeIndex |= 16;
-				if (CubeValues[5] < 0.0f) CubeIndex |= 32;
-				if (CubeValues[6] < 0.0f) CubeIndex |= 64;
-				if (CubeValues[7] < 0.0f) CubeIndex |= 128;
-
-				//Find the vertices where the surface intersects the cube
-				if (edgeTable[CubeIndex] & 1){
-					P = VertexInterp(C0,C1,CubeValues[0],CubeValues[1]);
-					VertexList[0] = P;
-					Q = C0;
-				}
-				if (edgeTable[CubeIndex] & 2){
-					P = VertexInterp(C1,C2,CubeValues[1],CubeValues[2]);
-					VertexList[1] = P;
-					Q = C1;
-				}
-				if (edgeTable[CubeIndex] & 4){
-					P = VertexInterp(C2,C3,CubeValues[2],CubeValues[3]);
-					VertexList[2] =	P;
-					Q = C2;
-				}
-				if (edgeTable[CubeIndex] & 8){
-					P = VertexInterp(C3,C0,CubeValues[3],CubeValues[0]);
-					VertexList[3] =	P;
-					Q = C3;
-				}
-				if (edgeTable[CubeIndex] & 16){
-					P = VertexInterp(C4,C5,CubeValues[4],CubeValues[5]);
-					VertexList[4] =	P;
-					Q = C4;
-				}
-				if (edgeTable[CubeIndex] & 32){
-					P = VertexInterp(C5,C6,CubeValues[5],CubeValues[6]);
-					VertexList[5] = P;
-					Q = C5;
-				}
-				if (edgeTable[CubeIndex] & 64){ 
-					P = VertexInterp(C6,C7,CubeValues[6],CubeValues[7]);
-					VertexList[6] = P;
-					Q = C6;
-				}
-				if (edgeTable[CubeIndex] & 128){
-					P = VertexInterp(C7,C4,CubeValues[7],CubeValues[4]);
-					VertexList[7] =	P;
-					Q = C7;
-				}
-				if (edgeTable[CubeIndex] & 256){
-					P = VertexInterp(C0,C4,CubeValues[0],CubeValues[4]);
-					VertexList[8] =	P;
-					Q = C0;
-				}
-				if (edgeTable[CubeIndex] & 512){
-					P = VertexInterp(C1,C5,CubeValues[1],CubeValues[5]);
-					VertexList[9] =	P;
-					Q = C1;
-				}
-				if (edgeTable[CubeIndex] & 1024){
-					P = VertexInterp(C2,C6,CubeValues[2],CubeValues[6]);
-					VertexList[10] = P;
-					Q = C2;
-				}
-				if (edgeTable[CubeIndex] & 2048){
-					P = VertexInterp(C3,C7,CubeValues[3],CubeValues[7]);
-					VertexList[11] = P;
-					Q = C3;
-				}
-
-				VertexCount=0;
-				for (int idx=0;idx<12;idx++)
-					LocalRemap[idx] = -1;
-
-				for (int idx=0;triTable[CubeIndex][idx]!=-1;idx++)
-				{
-					if(LocalRemap[triTable[CubeIndex][idx]] == -1)
-					{
-						NewVertexList[VertexCount] = VertexList[triTable[CubeIndex][idx]];
-						LocalRemap[triTable[CubeIndex][idx]] = VertexCount;
-						VertexCount++;
-					}
-				}
-
-				for (int idx=0;idx<VertexCount;idx++) {
-					P = NewVertexList[idx];
-					//P.x  += i;
-					//P.y  += j;
-					//P.z  += k;
-					cellvertices[idx] = P;
-				}
-
-				TriangleCount = 0;
-				for (int idx=0;triTable[CubeIndex][idx]!=-1;idx+=3) {
-					Triangles[TriangleCount][0] = LocalRemap[triTable[CubeIndex][idx+0]];
-					Triangles[TriangleCount][1] = LocalRemap[triTable[CubeIndex][idx+1]];
-					Triangles[TriangleCount][2] = LocalRemap[triTable[CubeIndex][idx+2]];
-					TriangleCount++;
-				}
-				int nTris = TriangleCount;
-
-				// Now add the local values to the DECL data structure
-				HalfEdge.resize(nTris*3);
-				int idx_edge=0;
-				for (int idx=0; idx<TriangleCount; idx++){
-					int V1 = Triangles[idx][0];
-					int V2 = Triangles[idx][1];
-					int V3 = Triangles[idx][2];
-					// first edge: V1->V2 
-					HalfEdge[idx_edge][0] = V1;  // first vertex
-					HalfEdge[idx_edge][1] = V2;  // second vertex
-					HalfEdge[idx_edge][2] = idx; // triangle
-					HalfEdge[idx_edge][3] = -1;  // twin
-					HalfEdge[idx_edge][4] = idx_edge+2;  // previous edge
-					HalfEdge[idx_edge][5] = idx_edge+1;  // next edge
-					idx_edge++;
-					// second edge: V2->V3
-					HalfEdge[idx_edge][0] = V2;  // first vertex
-					HalfEdge[idx_edge][1] = V3;  // second vertex
-					HalfEdge[idx_edge][2] = idx; // triangle
-					HalfEdge[idx_edge][3] = -1;  // twin
-					HalfEdge[idx_edge][4] = idx_edge-1;  // previous edge
-					HalfEdge[idx_edge][5] = idx_edge+1;  // next edge
-					idx_edge++;
-					// third edge: V3->V1
-					HalfEdge[idx_edge][0] = V3;  // first vertex
-					HalfEdge[idx_edge][1] = V1;  // second vertex
-					HalfEdge[idx_edge][2] = idx; // triangle
-					HalfEdge[idx_edge][3] = -1;  // twin
-					HalfEdge[idx_edge][4] = idx_edge-1;  // previous edge
-					HalfEdge[idx_edge][5] = idx_edge-2;  // next edge
-					idx_edge++;
-				}
-				int EdgeCount=idx_edge;
-				for (int idx=0; idx<EdgeCount; idx++){
-					int V1=HalfEdge[idx][0];
-					int V2=HalfEdge[idx][1];
-					// Find all the twins within the cube
-					for (int jdx=0; idx<EdgeCount; jdx++){
-						if (HalfEdge[jdx][1] == V1 && HalfEdge[jdx][0] == V2){
-							// this is the pair
-							HalfEdge[idx][3] = jdx;
-							HalfEdge[jdx][3] = idx;
-						}
-						if (HalfEdge[jdx][1] == V2 && HalfEdge[jdx][0] == V1 && !(idx==jdx)){
-							std::printf("WARNING: half edges with identical orientation! \n");
-						}
-					}	
-					// Use "ghost" twins if edge is on a cube face
-					P = cellvertices[V1];
-					Q = cellvertices[V2];
-					if (P.x == 0.0 && Q.x == 0.0) HalfEdge[idx_edge][3] = -1;  // ghost twin for x=0 face
-					if (P.x == 1.0 && Q.x == 1.0) HalfEdge[idx_edge][3] = -2;  // ghost twin for x=1 face
-					if (P.y == 0.0 && Q.y == 0.0) HalfEdge[idx_edge][3] = -3;  // ghost twin for y=0 face
-					if (P.y == 1.0 && Q.y == 1.0) HalfEdge[idx_edge][3] = -4;  // ghost twin for y=1 face
-					if (P.z == 0.0 && Q.z == 0.0) HalfEdge[idx_edge][3] = -5;  // ghost twin for z=0 face
-					if (P.z == 1.0 && Q.z == 1.0) HalfEdge[idx_edge][3] = -6;  // ghost twin for z=1 face
-				}
-				// Find all the angles
-				/*for (int idx=0; idx<EdgeCount; idx++){
-					int V1=HalfEdge[idx][0];
-					int V2=HalfEdge[idx][1];
-					int T1= HalfEdge[idx_edge][2];
-					int twin=HalfEdge[idx_edge][3];
-					if (twin == -1){
-						
-					}
-				}*/
-
-				// Map vertices to global coordinates
-				for (int idx=0;idx<VertexCount;idx++) {
-					P = cellvertices[idx];
-					P.x  += i;
-					P.y  += j;
-					P.z  += k;
-					cellvertices[idx] = P;
-				}
-			}
-		}
-	}
+  DECL object;
+  int e1,e2,e3;
+  FILE *TRIANGLES;
+  TRIANGLES = fopen("isosurface.stl","w");
+  fprintf(TRIANGLES,"solid isosurface\n");
+  int Nx = Field.size(0);
+  int Ny = Field.size(1);
+  int Nz = Field.size(2);
+  for (int k=1; k<Nz-1; k++){
+    for (int j=1; j<Ny-1; j++){
+      for (int i=1; i<Nx-1; i++){
+        object.LocalIsosurface(Field,isovalue,i,j,k);
+        for (int idx=0; idx<object.TriangleCount; idx++){
+          e1 = object.Face(idx);
+          e2 = object.halfedge.next(e1);
+          e3 = object.halfedge.next(e2);
+          auto P1 = object.vertex.coords(object.halfedge.v1(e1));
+          auto P2 = object.vertex.coords(object.halfedge.v1(e2));
+          auto P3 = object.vertex.coords(object.halfedge.v1(e3));
+          P1.x += 1.0*i; P1.y += 1.0*j; P1.z +=1.0*k;
+          P2.x += 1.0*i; P2.y += 1.0*j; P2.z +=1.0*k;
+          P3.x += 1.0*i; P3.y += 1.0*j; P3.z +=1.0*k;
+          fprintf(TRIANGLES,"vertex %f %f %f\n",P1.x,P1.y,P1.z);
+          fprintf(TRIANGLES,"vertex %f %f %f\n",P2.x,P2.y,P2.z);
+          fprintf(TRIANGLES,"vertex %f %f %f\n",P3.x,P3.y,P3.z);
 
+        }
+      }
+    }
+  }
+  fclose(TRIANGLES);
 }
-
-
-
diff --git a/analysis/dcel.h b/analysis/dcel.h
index c1c67178..0e6b2b1d 100644
--- a/analysis/dcel.h
+++ b/analysis/dcel.h
@@ -1,3 +1,6 @@
+#ifndef DCEL_INC
+#define DCEL_INC
+
 #include <vector>
 #include "analysis/pmmc.h"
 
@@ -79,3 +82,7 @@ public:
 private:
 	std::vector<int> FaceData;
 };
+
+void iso_surface(const Array<double>&Field, const double isovalue);
+
+#endif