Make the ADMMass computation working also with PUGH.

git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/ADMMass/trunk@20 54511f98-0e4f-0410-826e-eb8b393f5a1e

1 files changed, 54 insertions, 56 deletions

diff --git a/src/volume_integral.c b/src/volume_integral.c
index dcbf1a0..47374b3 100644
--- a/src/volume_integral.c
+++ b/src/volume_integral.c
@@ -194,29 +194,21 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
       }
 
 
-              /* Carpet does the following itself, but pugh does not. */
+    /* Carpet does the following itself, but pugh does not. */
 
     if (CCTK_IsThornActive("PUGH"))     
     {
-        /* Recompute the integrand on the symmetry and physical boundaries:
-           devide by 2 on the boundary faces cctk_lbnd values start from zero,
-           so must add one to get the total number of points.
-           If cctk_ubnd[2]=n, it means that it is the (n+1)-st point
-           (C notation) */
-
-        /* for short hand, right and left modified stencil values */ 
+      /* for short hand, right and left modified stencil values */ 
         const int lst  = ghost;
+        const int rst0 = cctk_lsh[0] - ghost;
         const int rst1 = cctk_lsh[1] - ghost;
         const int rst2 = cctk_lsh[2] - ghost;
-        const int rst3 = cctk_lsh[3] - ghost;
                 
         /*find the value of the "avoid_origin"  parameter */
         const void *avoid_origin_parameter_ptr =
-            CCTK_ParameterValString("avoid_origin","CartGrid3D");
+            CCTK_ParameterGet("avoid_origin","CartGrid3D",NULL);
         assert( avoid_origin_parameter_ptr != NULL );
         avoid_origin_parameter = *(const CCTK_INT *)avoid_origin_parameter_ptr;
-      
-        printf("avoid origin %d\n",avoid_origin_parameter);
     
         /* if  avoid_origin = yes (default), then do not devide by two. 
            This gives the correct result for the trapezoidal integration rule
@@ -226,51 +218,57 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
         
         if (! avoid_origin_parameter)
         {        
-            if (cctk_lbnd[2] == 0)
-                for(j = lst; j <= rst2; j++)
-                 for(i = lst; i <= rst1; i++)
-                 {
-                     ijk = CCTK_GFINDEX3D(cctkGH, i, j, lst);
-                     ADMMass_VolumeMass_GF[ijk] *= 0.5;
-                 }
-            if (cctk_ubnd[2]+1 == cctk_gsh[2])
-                for(j = lst; j <= rst2; j++)
-                 for(i = lst; i <= rst1; i++)
-                 {
-                     ijk = CCTK_GFINDEX3D(cctkGH, i, j, rst3);
-                     ADMMass_VolumeMass_GF[ijk] *= 0.5;
-                 }       
-            if (cctk_lbnd[1] == 0)
-                for(k = lst; k<= rst3; k++)
-                 for(i = lst; i<= rst1; i++)
-                 {
-                     ijk = CCTK_GFINDEX3D(cctkGH, i, lst, k);
-                     ADMMass_VolumeMass_GF[ijk] *= 0.5;
-                 }
-            if (cctk_ubnd[1]+1 == cctk_gsh[1])
-                for(k = lst; k<= rst3; k++)
-                 for(i = lst; i<= rst1; i++)
-                 {
-                     ijk = CCTK_GFINDEX3D(cctkGH, i, rst2, k);
-                     ADMMass_VolumeMass_GF[ijk] *= 0.5;
-                 }
-            if (cctk_lbnd[0] == 0)
-                for(k = lst; k <= rst3; k++)
-                 for(j = lst; j <= rst2; j++)
-                 {
-                     ijk = CCTK_GFINDEX3D(cctkGH, lst, j, k);
-                     ADMMass_VolumeMass_GF[ijk] *= 0.5;
-                 }
-            if (cctk_ubnd[0]+1 == cctk_gsh[0])
-                for(k = lst; k <= rst3; k++)
-                 for(j = lst; j <= rst2; j++)
-                 {
-                     ijk = CCTK_GFINDEX3D(cctkGH, rst1, j, k);
-                     ADMMass_VolumeMass_GF[ijk] *= 0.5;
-                 }
-          } /* end if avoid_origin */
+	  if (cctk_lbnd[2] == 0)
+	    for(j = lst; j <= rst1; j++)
+	      for(i = lst; i <= rst0; i++)
+		{
+		  ijk = CCTK_GFINDEX3D(cctkGH, i, j, lst);
+		  ADMMass_VolumeMass_GF[ijk] *= 0.5;
+		}
+	  /* Recompute the integrand on the symmetry and physical boundaries:
+	     devide by 2 on the boundary faces cctk_lbnd values start from zero,
+	     so must add one to get the total number of points.
+	     If cctk_ubnd[2]=n, it means that it is the (n+1)-st point
+	     (C notation) */
+	  
+	  if (cctk_ubnd[2]+1 == cctk_gsh[2])
+	    for(j = lst; j <= rst1; j++)
+	      for(i = lst; i <= rst0; i++)
+		{
+		  ijk = CCTK_GFINDEX3D(cctkGH, i, j, rst2);
+		  ADMMass_VolumeMass_GF[ijk] *= 0.5;
+		}       
+	  if (cctk_lbnd[1] == 0)
+	    for(k = lst; k<= rst2; k++)
+	      for(i = lst; i<= rst0; i++)
+		{
+		  ijk = CCTK_GFINDEX3D(cctkGH, i, lst, k);
+		  ADMMass_VolumeMass_GF[ijk] *= 0.5;
+		}
+	  if (cctk_ubnd[1]+1 == cctk_gsh[1])
+	    for(k = lst; k<= rst2; k++)
+	      for(i = lst; i<= rst0; i++)
+		{
+		  ijk = CCTK_GFINDEX3D(cctkGH, i, rst1, k);
+		  ADMMass_VolumeMass_GF[ijk] *= 0.5;
+		}
+	  if (cctk_lbnd[0] == 0)
+	    for(k = lst; k <= rst2; k++)
+	      for(j = lst; j <= rst1; j++)
+		{
+		  ijk = CCTK_GFINDEX3D(cctkGH, lst, j, k);
+		  ADMMass_VolumeMass_GF[ijk] *= 0.5;
+		}
+	  if (cctk_ubnd[0]+1 == cctk_gsh[0])
+	    for(k = lst; k <= rst2; k++)
+	      for(j = lst; j <= rst1; j++)
+		{
+		  ijk = CCTK_GFINDEX3D(cctkGH, rst0, j, k);
+		  ADMMass_VolumeMass_GF[ijk] *= 0.5;
+		}
+	} /* end if avoid_origin */
     } /* end if PUGH */
-
+    
     *grid_spacing_product = cctk_delta_space[0]*
                             cctk_delta_space[1]*
                             cctk_delta_space[2];