1 files changed, 31 insertions, 30 deletions

diff --git a/src/NuSurfacer.c b/src/NuSurfacer.c
index d12aba8..8e12186 100644
--- a/src/NuSurfacer.c
+++ b/src/NuSurfacer.c
@@ -3,8 +3,7 @@
 #endif
 #include <stdlib.h>
 #include <stdio.h>
-#include "common.h"
-#include "cctk.h"
+#include "IsoSurfacerGH.h"
 
 #include "CactusPUGH/PUGH/src/include/pugh.h"
 #include "NuSurfacer.h"
@@ -289,14 +288,24 @@ static TRIANGLE_CASES triCases[] = {
 		     for(nvals+=lst->nelem;lst->nelem<nvals;lst->nelem++) \
 	                lst->array[lst->nelem] = vals[lst->nelem]; }
 
-#define ComputePoints(pughGH, pindex,p1,p2) { int pointindex=pindex[0]; \
-p1[0]=(isofl)(i+xMapOffset[pointindex])*pughGH->dx0+pughGH->lx0;\
-p1[1]=(isofl)(j+yMapOffset[pointindex])*pughGH->dy0+pughGH->ly0;\
-p1[2]=(isofl)(k+zMapOffset[pointindex])*pughGH->dz0+pughGH->lz0;\
+#define ComputePoints(xcoords,ycoords,zcoords,pindex,p1,p2,xl,yl) { \
+int pointindex=pindex[0]; \
+int _ii,_jj,_kk,_idx;\
+_ii=i+xMapOffset[pointindex];\
+_jj=j+yMapOffset[pointindex];\
+_kk=k+zMapOffset[pointindex];\
+_idx=_ii+_jj*xl+_kk*yl;\
+p1[0]=(CCTK_REAL)xcoords[_idx];\
+p1[1]=(CCTK_REAL)ycoords[_idx];\
+p1[2]=(CCTK_REAL)zcoords[_idx];\
 pointindex=pindex[1];\
-p2[0]=(isofl)(i+xMapOffset[pointindex])*pughGH->dx0+pughGH->lx0;\
-p2[1]=(isofl)(j+yMapOffset[pointindex])*pughGH->dy0+pughGH->ly0;\
-p2[2]=(isofl)(k+zMapOffset[pointindex])*pughGH->dz0+pughGH->lz0;}
+_ii=i+xMapOffset[pointindex];\
+_jj=j+yMapOffset[pointindex];\
+_kk=k+zMapOffset[pointindex];\
+_idx=_ii+_jj*xl+_kk*yl;\
+p2[0]=(CCTK_REAL)xcoords[_idx];\
+p2[1]=(CCTK_REAL)ycoords[_idx];\
+p2[2]=(CCTK_REAL)zcoords[_idx];}
 
 #define SetEdgeMaskP(idx_,data_,value_,mask_) { \
 	register int _cmask=1; \
@@ -312,13 +321,13 @@ p2[2]=(isofl)(k+zMapOffset[pointindex])*pughGH->dz0+pughGH->lz0;}
 }
 
 #define AddInterpolatedVertex(isoval,idx,pindex,p1,p2,vertlist) { \
-  isofl _v,_d,_p1val,_p2val; int _i; \
+  CCTK_REAL4 _v,_d,_p1val,_p2val; int _i; \
   _p1val = data[idx+cellOffset[pindex[0]]]; \
   _p2val = data[idx+cellOffset[pindex[1]]]; \
   _d=(isoval-_p1val)/(_p2val-_p1val); \
   for( _i=0;_i<3;_i++){ \
      _v=(p2[_i]-p1[_i])*_d + p1[_i]; \
-     ArrayAddElem(isofl,vertlist,_v); \
+     ArrayAddElem(CCTK_REAL4,vertlist,_v); \
   } \
 }
 
@@ -332,7 +341,7 @@ typedef struct IsoCell {
 } IsoCell;
 
 IsoCell *cellslice;
-DeclareArray(isofl,FloatArray); /* declare a list of floats for vertices */
+DeclareArray(CCTK_REAL4,FloatArray); /* declare a list of floats for vertices */
 DeclareArray(CCTK_INT4,IntArray); /* declare a list of intes for triangle connectivity */
 FloatArray *vertexlist;
 IntArray *trianglelist;
@@ -361,14 +370,13 @@ static int isoCellVindex[12]={1,0,1,0,
 			      1,0,1,0,
 			      2,2,2,2};
 							
-void InitCellArrays(cGH *GH){
+void InitCellArrays(nx,ny){
   int i,j;
-  int nx=GH->cctk_lsh [0],ny=GH->cctk_lsh [1];
   int ix=1,iy=nx,iz=nx*ny;
   cellslice = New(IsoCell,2*nx*ny); /* 2 slices */
   for(i=0;i<nx*ny*2;i++) 
     for(j=0;j<3;j++) cellslice[i].v[j]=-1; /* init all verts */
-  vertexlist = New(FloatArray,1); InitArray(isofl,vertexlist);
+  vertexlist = New(FloatArray,1); InitArray(CCTK_REAL4,vertexlist);
   trianglelist = New(IntArray,1); InitArray(CCTK_INT4,trianglelist);
   /* celloffset */
   cellOffset[0] = 0;
@@ -388,31 +396,24 @@ void InitCellArrays(cGH *GH){
   isoCellOffset[10]=-(ix);   isoCellOffset[11]=0;/* -yx,-y,-x,0*/
 }
 
-void NuFindSurface(cGH *GH, int index,CCTK_REAL isovalue,polypatch *results){
+void NuFindSurface(CCTK_REAL *data, 
+		   int nx,int ny,int nz,
+		   CCTK_REAL *xcoords, CCTK_REAL *ycoords, CCTK_REAL *zcoords,
+		   CCTK_REAL isovalue, polypatch *results){
   static int initialized_nusurf=0;
   TRIANGLE_CASES *triCase;
   EDGE_LIST  *edge;
   int i,j,k,imax,jmax,kmax;
   int timelevel;
   int idx,nverts=0; /* vertex and triangle counters for offsets */
-  CCTK_REAL *data;
-  int nx=GH->cctk_lsh [0],ny=GH->cctk_lsh [1],nz=GH->cctk_lsh [2],slicesize,slicemodulo;
-  pGH *pughGH;
-
-  pughGH = PUGH_pGH (GH);
+  int slicesize,slicemodulo;
 
   if(!initialized_nusurf){
     initialized_nusurf=1;
-    InitCellArrays(GH);
+    InitCellArrays(nx,ny);
   }
   ResetArray(vertexlist);   /* zero the vertex count (but keep storage) */
   ResetArray(trianglelist); /* zero the triangle count (but keep storage) */
-
-  timelevel = CCTK_NumTimeLevelsFromVarI (index);
-  if (timelevel > 0)
-    timelevel--;
-  data = (CCTK_REAL *) GH->data [index][timelevel];
-
   /*****Sliceinit */
   slicesize=nx*ny; slicemodulo=slicesize*2;
   if(!cellslice) cellslice = New(IsoCell,slicemodulo);
@@ -426,7 +427,7 @@ void NuFindSurface(cGH *GH, int index,CCTK_REAL isovalue,polypatch *results){
     int isocellidx;
     int mask;
     int edgenumber;
-    isofl p1[3],p2[3]; 
+    CCTK_REAL4 p1[3],p2[3]; 
     int isocellnum;
     int isocellvindex; /* 12 maps to 3 */
     int elem,*elemp;
@@ -455,7 +456,7 @@ void NuFindSurface(cGH *GH, int index,CCTK_REAL isovalue,polypatch *results){
 	  elem=*elemp;
 	  if(elem<0){ /* add that vertex */
 	    elem=*elemp=nverts;
-	    ComputePoints(pughGH, edges[edgenumber],p1,p2);
+	    ComputePoints(xcoords,ycoords,zcoords, edges[edgenumber],p1,p2,nx,slicesize);
 	    AddInterpolatedVertex(isovalue,idx,edges[edgenumber],p1,p2,vertexlist);
 	    nverts++;
 	  }