Added gradient normal computations for the Renderer Thorn.

Also reenabled steering for isovalue because it works.


git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGHIO/IsoSurfacer/trunk@21 bfcf8e34-485d-4d46-a995-1fd6fa6fb178

5 files changed, 72 insertions, 5 deletions

diff --git a/src/IsoSurfacer.c b/src/IsoSurfacer.c
index 6c20659..c41c745 100644
--- a/src/IsoSurfacer.c
+++ b/src/IsoSurfacer.c
@@ -196,7 +196,7 @@ static void computeIso(int index, cGH *GH, isosurfacerGH *myGH)
   i = CCTK_CoordIndex (-1,"z","cart3d");
   zcoords = (CCTK_REAL *) CCTK_VarDataPtrI (GH, timelevel, i);
   /* Actually perform the Isosurfacing Operation */
-  NuFindSurface(data,nx,ny,nz,xcoords,ycoords,zcoords,myGH->isovalue,&(myGH->perprocessor));
+  NuFindSurface(data,nx,ny,nz,xcoords,ycoords,zcoords,myGH->isovalue,0,&(myGH->perprocessor));
   /* And collect the geometry to node 0 */
   CollectData(GH, &(myGH->perprocessor), &(myGH->totals));
   /* and for the writers, collect min,max to node 0 */
diff --git a/src/IsoSurfacerGH.h b/src/IsoSurfacerGH.h
index 9f68814..3791bfd 100644
--- a/src/IsoSurfacerGH.h
+++ b/src/IsoSurfacerGH.h
@@ -25,12 +25,15 @@ typedef union IsoCommand {
   CmdSet cmd;
 } IsoCommand;
 
+/* Should change name to IsoGeometry */
 typedef struct polypatch
 {
   CCTK_REAL4    *verts;
   CCTK_INT4      nverts;
   CCTK_INT4     *polys;
   CCTK_INT4      npolys;
+  CCTK_REAL4    *norms; /* these are vertex normals */
+  CCTK_INT4      nnorms;
 } polypatch;
 
 typedef struct isosurfacerGH
@@ -41,7 +44,8 @@ typedef struct isosurfacerGH
   short int   formats;
   short int   outfreq;
   short int   firstIteration;
-  short int   RunIsoSurfacer;
+  short int  RunIsoSurfacer;
+  short int  ComputeNormals;
   polypatch   perprocessor,totals; /* un-collected geometry for this node */
 } isosurfacerGH;
 
diff --git a/src/IsoSurfacerInit.c b/src/IsoSurfacerInit.c
index 36f2a0a..16c872d 100644
--- a/src/IsoSurfacerInit.c
+++ b/src/IsoSurfacerInit.c
@@ -24,6 +24,7 @@ void *IsoSurfacer_SetupGH (tFleshConfig *config,
   myGH->formats=0;
   myGH->outfreq=0;
   myGH->firstIteration=0;
+  myGH->ComputeNormals=0;
   myGH->isovalue=0.0;
   myGH->perprocessor.verts  = myGH->totals.verts  = NULL;
   myGH->perprocessor.nverts = myGH->totals.nverts = 0;
@@ -56,6 +57,7 @@ int IsoSurfacer_InitGH (cGH *GH){
   myGH->formats=0;
   myGH->outfreq=output_frequency;
   myGH->firstIteration=output_start;
+  myGH->ComputeNormals=compute_normals;
   myGH->isovalue=isovalue;
   myGH->perprocessor.verts  = myGH->totals.verts  = NULL;
   myGH->perprocessor.nverts = myGH->totals.nverts = 0;
diff --git a/src/NuSurfacer.c b/src/NuSurfacer.c
index 8e12186..c38c805 100644
--- a/src/NuSurfacer.c
+++ b/src/NuSurfacer.c
@@ -3,6 +3,7 @@
 #endif
 #include <stdlib.h>
 #include <stdio.h>
+#include <math.h>
 #include "IsoSurfacerGH.h"
 
 #include "CactusPUGH/PUGH/src/include/pugh.h"
@@ -320,6 +321,51 @@ p2[2]=(CCTK_REAL)zcoords[_idx];}
 	if(data_[idx_+cellOffset[7]]>=value_) mask_|=_cmask; \
 }
 
+/*
+  Need to get celloffset in the 3rd dimension.
+  Then normalize the gradient magnitude to 1.
+  (see if VTK has a better grad solution)
+
+  Need to determine coorect dim orientation.
+  pindex=?
+  celloffset=?
+  idx=?
+
+  cell vector orientation?  (or-together pindex and use to compute
+  from table the root node required to look in each vector direction)
+  well really, we can just take pindex[0] and look 3 directions and
+  then pindex[1] and look 3 directions.  Add the vectors weighted by
+  the distance between the two points then normalize.
+  Even dirtier, take the gradient at the lowest order point?
+
+  Either way, we need to bounds-check on this.
+  (urk... celloffset is insufficient)
+  We need to do celloff for each direction
+ */
+#define AddGradientNormal(isoval,idx,pindex,p1,p2,normlist) { \
+  CCTK_REAL4 _d,_p1val,_p2val,_p1vec[3],_p2vec[3],_scale; \
+  int _i,_celloff[3];\
+  _celloff[0]=cellOffset[1]; _celloff[1]=cellOffset[3]; _celloff[2]=cellOffset[4];\
+  _p1vec[0]=_p1vec[1]=_p1vec[2]=-data[idx+cellOffset[pindex[0]]];\
+  _p2vec[0]=_p2vec[1]=_p2vec[2]=-data[idx+cellOffset[pindex[1]]];\
+  _p1val = data[idx+cellOffset[pindex[0]]]; \
+  _p2val = data[idx+cellOffset[pindex[1]]]; \
+  _d=(isoval-_p1val)/(_p2val-_p1val); \
+ /* printf("_d=%g\n",_d); */\
+  for( _scale=0,_i=0;_i<3;_i++){ \
+     _p1vec[_i]+=data[idx+cellOffset[pindex[0]]+_celloff[_i]]; \
+     _p2vec[_i]+=data[idx+cellOffset[pindex[1]]+_celloff[_i]]; \
+     _p1vec[_i]=(_p2vec[_i]-_p1vec[_i])*_d + _p1vec[_i]; \
+     _scale+=(_p1vec[_i]*_p1vec[_i]); \
+  } \
+ /* printf("Accumulated _scale=%g\n",_scale);*/\
+  _scale=1.0/sqrt(_scale); \
+  for(_i=0;_i<3;_i++){ \
+      _p1vec[_i]*=_scale; \
+      ArrayAddElem(CCTK_REAL4,normlist,_p1vec[_i]); \
+  } \
+}
+
 #define AddInterpolatedVertex(isoval,idx,pindex,p1,p2,vertlist) { \
   CCTK_REAL4 _v,_d,_p1val,_p2val; int _i; \
   _p1val = data[idx+cellOffset[pindex[0]]]; \
@@ -343,7 +389,7 @@ typedef struct IsoCell {
 IsoCell *cellslice;
 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;
+FloatArray *vertexlist,*normallist;
 IntArray *trianglelist;
 static int edges[12][2] = { {0,1}, {1,2}, {3,2}, {0,3},
                             {4,5}, {5,6}, {7,6}, {4,7},
@@ -369,6 +415,8 @@ static int isoCellOffset[12]={0,0,0,0, /* -zy,-z,-z,-zx, */
 static int isoCellVindex[12]={1,0,1,0,
 			      1,0,1,0,
 			      2,2,2,2};
+
+/* cell vector orientation? */
 							
 void InitCellArrays(nx,ny){
   int i,j;
@@ -377,6 +425,7 @@ void InitCellArrays(nx,ny){
   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(CCTK_REAL4,vertexlist);
+  normallist = New(FloatArray,1); InitArray(CCTK_REAL4,normallist);
   trianglelist = New(IntArray,1); InitArray(CCTK_INT4,trianglelist);
   /* celloffset */
   cellOffset[0] = 0;
@@ -399,7 +448,7 @@ void InitCellArrays(nx,ny){
 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){
+		   CCTK_REAL isovalue, int compute_normals, polypatch *results){
   static int initialized_nusurf=0;
   TRIANGLE_CASES *triCase;
   EDGE_LIST  *edge;
@@ -413,6 +462,8 @@ void NuFindSurface(CCTK_REAL *data,
     InitCellArrays(nx,ny);
   }
   ResetArray(vertexlist);   /* zero the vertex count (but keep storage) */
+  if(compute_normals)
+    ResetArray(normallist);   /* zero the norm count (but keep storage) */
   ResetArray(trianglelist); /* zero the triangle count (but keep storage) */
   /*****Sliceinit */
   slicesize=nx*ny; slicemodulo=slicesize*2;
@@ -458,6 +509,8 @@ void NuFindSurface(CCTK_REAL *data,
 	    elem=*elemp=nverts;
 	    ComputePoints(xcoords,ycoords,zcoords, edges[edgenumber],p1,p2,nx,slicesize);
 	    AddInterpolatedVertex(isovalue,idx,edges[edgenumber],p1,p2,vertexlist);
+	    if(compute_normals)
+	      AddGradientNormal(isovalue,idx,edges[edgenumber],p1,p2,normallist);
 	    nverts++;
 	  }
 	  ArrayAddElem(CCTK_INT4,trianglelist,elem); /* this is just one vertex! */
@@ -469,5 +522,13 @@ void NuFindSurface(CCTK_REAL *data,
   results->verts = vertexlist->array;
   results->npolys = (trianglelist->nelem)?(trianglelist->nelem/3):(trianglelist->nelem);
   results->polys = trianglelist->array;
+  if(compute_normals){
+    results->nnorms = nverts; /* it is vertex-normals here */
+    results->norms = normallist->array;
+  }
+  else {
+    results->nnorms=0;
+    results->norms=NULL;
+  }
 }	 
   
diff --git a/src/NuSurfacer.h b/src/NuSurfacer.h
index f5c5c29..20d2cf5 100644
--- a/src/NuSurfacer.h
+++ b/src/NuSurfacer.h
@@ -5,6 +5,6 @@
 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);
+		   CCTK_REAL isovalue, int compute_normals, polypatch *results);
 
 #endif