The carpet patch timed Sat Mar 4 18:38:10 CET 2006, not allowing any longer the use of CCTK_DELTA_SPACE in global mode, had broken our way of computing volume integrals. Here is the upgrade.

Further changes are the transformation of divisions into moltiplications and the boundary corrections in the computation of the integral when PUGH is used. This last part is not tested (but will be), anyway who uses PUGH nowadays?


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

2 files changed, 177 insertions, 106 deletions

diff --git a/src/surface_integral.c b/src/surface_integral.c
index fad504b..b8a3998 100644
--- a/src/surface_integral.c
+++ b/src/surface_integral.c
@@ -43,9 +43,9 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
     DECLARE_CCTK_ARGUMENTS
     DECLARE_CCTK_PARAMETERS
 
-    CCTK_INT i,j,k, ijk, ierr, ghost;
+    CCTK_INT i,j,k, ijk, ierr, ghost, ti, tj, tk, tl;
     CCTK_INT  x_min_i, x_max_i, y_min_j, y_max_j, z_min_k, z_max_k;
-    CCTK_REAL detg, idetg, oneDX, oneDY, oneDZ, twoDX, twoDY, twoDZ, ds[3];
+    CCTK_REAL detg, idetg, ds[3];
     CCTK_REAL u[3][3], dg[3][3][3];
 
     CCTK_REAL physical_min[3];
@@ -58,18 +58,22 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
 
     char coordinate_parameter_type[9];
 
+#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h"
+
     /* grid-function strides for ADMMacros */
     const CCTK_INT di = 1;
     const CCTK_INT dj = cctk_lsh[0];
     const CCTK_INT dk = cctk_lsh[0]*cctk_lsh[1];
 
     /* deonminators for derivatives */
-    twoDX = 2.0 * CCTK_DELTA_SPACE(0);
-    twoDY = 2.0 * CCTK_DELTA_SPACE(1);
-    twoDZ = 2.0 * CCTK_DELTA_SPACE(2);
-    oneDX =  CCTK_DELTA_SPACE(0);
-    oneDY =  CCTK_DELTA_SPACE(1);
-    oneDZ =  CCTK_DELTA_SPACE(2);
+    const CCTK_REAL OneOverTwoDX = 1.0 / (2.0 * CCTK_DELTA_SPACE(0));
+    const CCTK_REAL OneOverTwoDY = 1.0 / (2.0 * CCTK_DELTA_SPACE(1));
+    const CCTK_REAL OneOverTwoDZ = 1.0 / (2.0 * CCTK_DELTA_SPACE(2));
+
+    const CCTK_REAL oneDX =  CCTK_DELTA_SPACE(0);
+    const CCTK_REAL oneDY =  CCTK_DELTA_SPACE(1);
+    const CCTK_REAL oneDZ =  CCTK_DELTA_SPACE(2);
+
 
     x_min_i = -INT_MAX;
     x_max_i =  INT_MAX; 
@@ -156,8 +160,6 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
     if ( (z[0] < *ADMMass_box_z_max) && (*ADMMass_box_z_max < z[ijk]) )
       z_max_k = find_closest(cctkGH, cctk_lsh, cctk_delta_space, ghost, z, *ADMMass_box_z_max, 2);
 
-#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h"
-
     for(i=ghost; i<cctk_lsh[0]-ghost; i++)
      for(j=ghost; j<cctk_lsh[1]-ghost; j++)
       for(k=ghost; k<cctk_lsh[2]-ghost; k++)
@@ -216,29 +218,29 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
                   u[2][1] = UPPERMET_UYZ;
                   u[2][2] = UPPERMET_UZZ;
 
-		  dg[0][0][0] = ( gxx[di+ijk] - gxx[-di+ijk] ) / twoDX;
-		  dg[0][0][1] = ( gxx[dj+ijk] - gxx[-dj+ijk] ) / twoDY;
-		  dg[0][0][2] = ( gxx[dk+ijk] - gxx[-dk+ijk] ) / twoDZ;
+		  dg[0][0][0] = ( gxx[di+ijk] - gxx[-di+ijk] ) * OneOverTwoDX;
+		  dg[0][0][1] = ( gxx[dj+ijk] - gxx[-dj+ijk] ) * OneOverTwoDY;
+		  dg[0][0][2] = ( gxx[dk+ijk] - gxx[-dk+ijk] ) * OneOverTwoDZ;
 		  
-		  dg[0][1][0] = ( gxy[di+ijk] - gxy[-di+ijk] ) / twoDX;
-		  dg[0][1][1] = ( gxy[dj+ijk] - gxy[-dj+ijk] ) / twoDY;
-		  dg[0][1][2] = ( gxy[dk+ijk] - gxy[-dk+ijk] ) / twoDZ;
+		  dg[0][1][0] = ( gxy[di+ijk] - gxy[-di+ijk] ) * OneOverTwoDX;
+		  dg[0][1][1] = ( gxy[dj+ijk] - gxy[-dj+ijk] ) * OneOverTwoDY;
+		  dg[0][1][2] = ( gxy[dk+ijk] - gxy[-dk+ijk] ) * OneOverTwoDZ;
 		  
-		  dg[0][2][0] = ( gxz[di+ijk] - gxz[-di+ijk] ) / twoDX;
-		  dg[0][2][1] = ( gxz[dj+ijk] - gxz[-dj+ijk] ) / twoDY;
-		  dg[0][2][2] = ( gxz[dk+ijk] - gxz[-dk+ijk] ) / twoDZ;
+		  dg[0][2][0] = ( gxz[di+ijk] - gxz[-di+ijk] ) * OneOverTwoDX;
+		  dg[0][2][1] = ( gxz[dj+ijk] - gxz[-dj+ijk] ) * OneOverTwoDY;
+		  dg[0][2][2] = ( gxz[dk+ijk] - gxz[-dk+ijk] ) * OneOverTwoDZ;
 		  
 		  dg[1][0][0] = dg[0][1][0];
 		  dg[1][0][1] = dg[0][1][1];
 		  dg[1][0][2] = dg[0][1][2];
 		  
-		  dg[1][1][0] = ( gyy[di+ijk] - gyy[-di+ijk] ) / twoDX;
-		  dg[1][1][1] = ( gyy[dj+ijk] - gyy[-dj+ijk] ) / twoDY;
-		  dg[1][1][2] = ( gyy[dk+ijk] - gyy[-dk+ijk] ) / twoDZ;
+		  dg[1][1][0] = ( gyy[di+ijk] - gyy[-di+ijk] ) * OneOverTwoDX;
+		  dg[1][1][1] = ( gyy[dj+ijk] - gyy[-dj+ijk] ) * OneOverTwoDY;
+		  dg[1][1][2] = ( gyy[dk+ijk] - gyy[-dk+ijk] ) * OneOverTwoDZ;
 		  
-		  dg[1][2][0] = ( gyz[di+ijk] - gyz[-di+ijk] ) / twoDX;
-		  dg[1][2][1] = ( gyz[dj+ijk] - gyz[-dj+ijk] ) / twoDY;
-		  dg[1][2][2] = ( gyz[dk+ijk] - gyz[-dk+ijk] ) / twoDZ;
+		  dg[1][2][0] = ( gyz[di+ijk] - gyz[-di+ijk] ) * OneOverTwoDX;
+		  dg[1][2][1] = ( gyz[dj+ijk] - gyz[-dj+ijk] ) * OneOverTwoDY;
+		  dg[1][2][2] = ( gyz[dk+ijk] - gyz[-dk+ijk] ) * OneOverTwoDZ;
 		  
 		  dg[2][0][0] = dg[0][2][0];
 		  dg[2][0][1] = dg[0][2][1];
@@ -248,14 +250,14 @@ void ADMMass_Surface(CCTK_ARGUMENTS)
 		  dg[2][1][1] = dg[1][2][1];
 		  dg[2][1][2] = dg[1][2][2];
 		  
-		  dg[2][2][0] = ( gzz[di+ijk] - gzz[-di+ijk] ) / twoDX;
-		  dg[2][2][1] = ( gzz[dj+ijk] - gzz[-dj+ijk] ) / twoDY;
-		  dg[2][2][2] = ( gzz[dk+ijk] - gzz[-dk+ijk] ) / twoDZ;
-
-                  for (int ti = 0; ti < 3; ti++)
-                   for (int tj = 0; tj < 3; tj++)
-                    for (int tk = 0; tk < 3; tk++)
-                     for (int tl = 0; tl < 3; tl++)
+		  dg[2][2][0] = ( gzz[di+ijk] - gzz[-di+ijk] ) * OneOverTwoDX;
+		  dg[2][2][1] = ( gzz[dj+ijk] - gzz[-dj+ijk] ) * OneOverTwoDY;
+		  dg[2][2][2] = ( gzz[dk+ijk] - gzz[-dk+ijk] ) * OneOverTwoDZ;
+
+                  for (ti = 0; ti < 3; ti++)
+                   for (tj = 0; tj < 3; tj++)
+                    for (tk = 0; tk < 3; tk++)
+                     for (tl = 0; tl < 3; tl++)
 		       {
                          ADMMass_SurfaceMass_GF[ijk] +=
 			   u[ti][tj] * u[tk][tl] * ( dg[ti][tk][tj] - dg[ti][tj][tk] ) * ds[tk];
@@ -288,7 +290,7 @@ void ADMMass_Surface_Global(CCTK_ARGUMENTS)
     ADMMass_SurfaceMass[*ADMMass_LoopCounter] /=  16.0*PI;
 
     CCTK_VInfo(CCTK_THORNSTRING,
-	       "detector %d: ADM mass  %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g\n",
+	       "detector %d: ADM mass  %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g",
 	       *ADMMass_LoopCounter,
 	       ADMMass_SurfaceMass[*ADMMass_LoopCounter],
 	       *ADMMass_box_x_min,*ADMMass_box_x_max,
@@ -345,7 +347,7 @@ void ADMMass_Surface_Lapse_Global(CCTK_ARGUMENTS)
     ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter] /=  16.0*PI;    
 
     CCTK_VInfo(CCTK_THORNSTRING,
-	       "detector %d: ADM mass with lapse  %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g\n",
+	       "detector %d: ADM mass with lapse  %g, surface (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g",
 	       *ADMMass_LoopCounter,
 	       ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter],
 	       *ADMMass_box_x_min,*ADMMass_box_x_max,
diff --git a/src/volume_integral.c b/src/volume_integral.c
index ee4009d..b64feaa 100644
--- a/src/volume_integral.c
+++ b/src/volume_integral.c
@@ -25,15 +25,16 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
     DECLARE_CCTK_ARGUMENTS
     DECLARE_CCTK_PARAMETERS
 
-    CCTK_INT i,j,k, ijk, ghost;
+    CCTK_INT i,j,k, ijk, ghost, ti, tj, tk, tl;
     CCTK_REAL detg, idetg;
     CCTK_REAL u[3][3], dg[3][3][3];
-    CCTK_REAL radius, twoDX, twoDY, twoDZ;
+    CCTK_REAL radius;
 
     CCTK_INT mask_descriptor, state_descriptor_outside;
 
-    mask_descriptor = SpaceMask_GetTypeBits("OutsideMask");
-    state_descriptor_outside = SpaceMask_GetStateBits("OutsideMask", "outside");
+    int avoid_origin_parameter;
+
+#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h"
 
     /* grid-function strides */
     const CCTK_INT di = 1;
@@ -41,9 +42,12 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
     const CCTK_INT dk = cctk_lsh[0]*cctk_lsh[1];
 
     /* deonminators for derivatives */
-    twoDX = 2.0 * CCTK_DELTA_SPACE(0);
-    twoDY = 2.0 * CCTK_DELTA_SPACE(1);
-    twoDZ = 2.0 * CCTK_DELTA_SPACE(2);
+    const CCTK_REAL OneOverTwoDX = 1.0 / (2.0 * CCTK_DELTA_SPACE(0));
+    const CCTK_REAL OneOverTwoDY = 1.0 / (2.0 * CCTK_DELTA_SPACE(1));
+    const CCTK_REAL OneOverTwoDZ = 1.0 / (2.0 * CCTK_DELTA_SPACE(2));
+
+    mask_descriptor = SpaceMask_GetTypeBits("OutsideMask");
+    state_descriptor_outside = SpaceMask_GetStateBits("OutsideMask", "outside");
 
     if (ADMMass_use_surface_distance_as_volume_radius &&
         (ADMMass_volume_radius[*ADMMass_LoopCounter] < 0.0))
@@ -71,8 +75,6 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
 	ghost = *(const int *)ghost_ptr;
       }
 
-#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h"
-
     for(i=1; i<cctk_lsh[0]-1; i++)
      for(j=1; j<cctk_lsh[1]-1; j++)
       for(k=1; k<cctk_lsh[2]-1; k++)   
@@ -94,29 +96,29 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
           u[2][1] = UPPERMET_UYZ;
           u[2][2] = UPPERMET_UZZ;
 
-	  dg[0][0][0] = ( gxx[di+ijk] - gxx[-di+ijk] ) / twoDX;
-          dg[0][0][1] = ( gxx[dj+ijk] - gxx[-dj+ijk] ) / twoDY;
-          dg[0][0][2] = ( gxx[dk+ijk] - gxx[-dk+ijk] ) / twoDZ;
+	  dg[0][0][0] = ( gxx[di+ijk] - gxx[-di+ijk] ) * OneOverTwoDX;
+          dg[0][0][1] = ( gxx[dj+ijk] - gxx[-dj+ijk] ) * OneOverTwoDY;
+          dg[0][0][2] = ( gxx[dk+ijk] - gxx[-dk+ijk] ) * OneOverTwoDZ;
 
-          dg[0][1][0] = ( gxy[di+ijk] - gxy[-di+ijk] ) / twoDX;
-          dg[0][1][1] = ( gxy[dj+ijk] - gxy[-dj+ijk] ) / twoDY;
-          dg[0][1][2] = ( gxy[dk+ijk] - gxy[-dk+ijk] ) / twoDZ;
+          dg[0][1][0] = ( gxy[di+ijk] - gxy[-di+ijk] ) * OneOverTwoDX;
+          dg[0][1][1] = ( gxy[dj+ijk] - gxy[-dj+ijk] ) * OneOverTwoDY;
+          dg[0][1][2] = ( gxy[dk+ijk] - gxy[-dk+ijk] ) * OneOverTwoDZ;
 
-          dg[0][2][0] = ( gxz[di+ijk] - gxz[-di+ijk] ) / twoDX;
-          dg[0][2][1] = ( gxz[dj+ijk] - gxz[-dj+ijk] ) / twoDY;
-          dg[0][2][2] = ( gxz[dk+ijk] - gxz[-dk+ijk] ) / twoDZ;
+          dg[0][2][0] = ( gxz[di+ijk] - gxz[-di+ijk] ) * OneOverTwoDX;
+          dg[0][2][1] = ( gxz[dj+ijk] - gxz[-dj+ijk] ) * OneOverTwoDY;
+          dg[0][2][2] = ( gxz[dk+ijk] - gxz[-dk+ijk] ) * OneOverTwoDZ;
 
           dg[1][0][0] = dg[0][1][0];
           dg[1][0][1] = dg[0][1][1];
           dg[1][0][2] = dg[0][1][2];
 
-          dg[1][1][0] = ( gyy[di+ijk] - gyy[-di+ijk] ) / twoDX;
-          dg[1][1][1] = ( gyy[dj+ijk] - gyy[-dj+ijk] ) / twoDY;
-          dg[1][1][2] = ( gyy[dk+ijk] - gyy[-dk+ijk] ) / twoDZ;
+          dg[1][1][0] = ( gyy[di+ijk] - gyy[-di+ijk] ) * OneOverTwoDX;
+          dg[1][1][1] = ( gyy[dj+ijk] - gyy[-dj+ijk] ) * OneOverTwoDY;
+          dg[1][1][2] = ( gyy[dk+ijk] - gyy[-dk+ijk] ) * OneOverTwoDZ;
 
-          dg[1][2][0] = ( gyz[di+ijk] - gyz[-di+ijk] ) / twoDX;
-          dg[1][2][1] = ( gyz[dj+ijk] - gyz[-dj+ijk] ) / twoDY;
-          dg[1][2][2] = ( gyz[dk+ijk] - gyz[-dk+ijk] ) / twoDZ;
+          dg[1][2][0] = ( gyz[di+ijk] - gyz[-di+ijk] ) * OneOverTwoDX;
+          dg[1][2][1] = ( gyz[dj+ijk] - gyz[-dj+ijk] ) * OneOverTwoDY;
+          dg[1][2][2] = ( gyz[dk+ijk] - gyz[-dk+ijk] ) * OneOverTwoDZ;
 
           dg[2][0][0] = dg[0][2][0];
           dg[2][0][1] = dg[0][2][1];
@@ -125,42 +127,14 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
           dg[2][1][0] = dg[1][2][0];
           dg[2][1][1] = dg[1][2][1];
           dg[2][1][2] = dg[1][2][2];
-
-          dg[2][2][0] = ( gzz[di+ijk] - gzz[-di+ijk] ) / twoDX;
-          dg[2][2][1] = ( gzz[dj+ijk] - gzz[-dj+ijk] ) / twoDY;
-          dg[2][2][2] = ( gzz[dk+ijk] - gzz[-dk+ijk] ) / twoDZ;
-
-	  /*	  dg[0][0][0] = DXDG_DXDGXX;
-          dg[0][0][1] = DYDG_DYDGXX;
-          dg[0][0][2] = DZDG_DZDGXX;
-          dg[0][1][0] = DXDG_DXDGXY;
-          dg[0][1][1] = DYDG_DYDGXY;
-          dg[0][1][2] = DZDG_DZDGXY;
-          dg[0][2][0] = DXDG_DXDGXZ;
-          dg[0][2][1] = DYDG_DYDGXZ;
-          dg[0][2][2] = DZDG_DZDGXZ;
-          dg[1][0][0] = DXDG_DXDGXY;
-          dg[1][0][1] = DYDG_DYDGXY;
-          dg[1][0][2] = DZDG_DZDGXY;
-          dg[1][1][0] = DXDG_DXDGYY;
-          dg[1][1][1] = DYDG_DYDGYY;
-          dg[1][1][2] = DZDG_DZDGYY;
-          dg[1][2][0] = DXDG_DXDGYZ;
-          dg[1][2][1] = DYDG_DYDGYZ;
-          dg[1][2][2] = DZDG_DZDGYZ;
-          dg[2][0][0] = DXDG_DXDGXZ;
-          dg[2][0][1] = DYDG_DYDGXZ;
-          dg[2][0][2] = DZDG_DZDGXZ;
-          dg[2][1][0] = DXDG_DXDGYZ;
-          dg[2][1][1] = DYDG_DYDGYZ;
-          dg[2][1][2] = DZDG_DZDGYZ;
-          dg[2][2][0] = DXDG_DXDGZZ;
-          dg[2][2][1] = DYDG_DYDGZZ;
-          dg[2][2][2] = DZDG_DZDGZZ;*/
-
-          for (int ti = 0; ti < 3; ti++)
-           for (int tj = 0; tj < 3; tj++)
-            for (int tk = 0; tk < 3; tk++)
+	  
+          dg[2][2][0] = ( gzz[di+ijk] - gzz[-di+ijk] ) * OneOverTwoDX;
+          dg[2][2][1] = ( gzz[dj+ijk] - gzz[-dj+ijk] ) * OneOverTwoDY;
+          dg[2][2][2] = ( gzz[dk+ijk] - gzz[-dk+ijk] ) * OneOverTwoDZ;
+
+          for (ti = 0; ti < 3; ti++)
+           for (tj = 0; tj < 3; tj++)
+            for (tk = 0; tk < 3; tk++)
             {
                 ADMMass_VolumeMass_pot_x[ijk] +=
                     u[ti][tj] * u[tk][0] *
@@ -198,16 +172,111 @@ void ADMMass_Volume(CCTK_ARGUMENTS)
           {
               ADMMass_VolumeMass_GF[ijk] =
                ((ADMMass_VolumeMass_pot_x[CCTK_GFINDEX3D(cctkGH,i+1,j,k)]-
-                      ADMMass_VolumeMass_pot_x[CCTK_GFINDEX3D(cctkGH,i-1,j,k)])/twoDX+
+                      ADMMass_VolumeMass_pot_x[CCTK_GFINDEX3D(cctkGH,i-1,j,k)])*OneOverTwoDX+
                 (ADMMass_VolumeMass_pot_y[CCTK_GFINDEX3D(cctkGH,i,j+1,k)]-
-                      ADMMass_VolumeMass_pot_y[CCTK_GFINDEX3D(cctkGH,i,j-1,k)])/twoDY+
+                      ADMMass_VolumeMass_pot_y[CCTK_GFINDEX3D(cctkGH,i,j-1,k)])*OneOverTwoDY+
                 (ADMMass_VolumeMass_pot_z[CCTK_GFINDEX3D(cctkGH,i,j,k+1)]-
-                      ADMMass_VolumeMass_pot_z[CCTK_GFINDEX3D(cctkGH,i,j,k-1)])/twoDZ);
+                      ADMMass_VolumeMass_pot_z[CCTK_GFINDEX3D(cctkGH,i,j,k-1)])*OneOverTwoDZ);
+	  }
+      }
+
 
+	      /* 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 */ 
+	const int lst  = 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");
+	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 on the symmetry boundaries, 
+	   while (even if it is wrong not to devide by two the points on the faces that are 
+	   physical outer boundaries) it should not spoil the result. */
+	
+	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] = 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] = 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] = 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] = 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] = 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] = 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];
+     
 #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_undefine.h"
-}
+}      
 
 void ADMMass_Volume_Global(CCTK_ARGUMENTS)
 {
@@ -226,14 +295,14 @@ void ADMMass_Volume_Global(CCTK_ARGUMENTS)
                     &ADMMass_VolumeMass[*ADMMass_LoopCounter], 1,
                     CCTK_VarIndex("ADMMass::ADMMass_VolumeMass_GF")))
         CCTK_WARN(0, "Error while reducing ADMMass_VolumeMass_GF");
-    
+
     ADMMass_VolumeMass[*ADMMass_LoopCounter] *=
-      cctk_delta_space[0]*cctk_delta_space[1]*cctk_delta_space[2] / (16.0*PI);
+      *grid_spacing_product / (16.0*PI);
 
     if (ADMMass_use_surface_distance_as_volume_radius)
       {
 	CCTK_VInfo(CCTK_THORNSTRING,
-		   "detector %d: ADM mass  %g, volume (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g\n",
+		   "detector %d: ADM mass  %g, volume (cube): xmin %g, xmax %g, ymin %g, ymax %g, zmin %g, zmax %g",
 		   *ADMMass_LoopCounter,
 		   ADMMass_SurfaceMass[*ADMMass_LoopCounter],
 		   *ADMMass_box_x_min,*ADMMass_box_x_max,
@@ -242,14 +311,14 @@ void ADMMass_Volume_Global(CCTK_ARGUMENTS)
       }
     else if (ADMMass_use_all_volume_as_volume_radius)
       {
-	    CCTK_VInfo(CCTK_THORNSTRING," detector %d: ADM mass %g, volume: the whole grid\n", 
+	    CCTK_VInfo(CCTK_THORNSTRING," detector %d: ADM mass %g, volume: the whole grid",
 	       *ADMMass_LoopCounter, ADMMass_VolumeMass[*ADMMass_LoopCounter]);
       }
     else
       {
-	CCTK_VInfo(CCTK_THORNSTRING," detector %d: ADM mass %g, volume (sphere of radius) %g\n", 
+	CCTK_VInfo(CCTK_THORNSTRING," detector %d: ADM mass %g, volume (sphere of radius) %g",
 	       *ADMMass_LoopCounter, ADMMass_VolumeMass[*ADMMass_LoopCounter],
 	       ADMMass_volume_radius[*ADMMass_LoopCounter]);
       }
-}
 
+}