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

4 files changed, 188 insertions, 107 deletions

diff --git a/interface.ccl b/interface.ccl
index c71e1f8..384c54d 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -47,3 +47,5 @@ CCTK_REAL  ADMMass_box type = scalar tags='checkpoint="no"'
  ADMMass_box_z_min
  ADMMass_box_z_max
 } "Physical coordinates of the surface on which the integral is computed"
+
+real grid_spacing_product type=SCALAR tags='checkpoint="no"' "product of cctk_delta_space, to be computed in local mode and later used in global mode (carpet problems)"
diff --git a/schedule.ccl b/schedule.ccl
index 42ebdd0..fc54cb1 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -22,11 +22,19 @@ schedule ADMMass_Loop IN ADMMass BEFORE ADMMass_Local
   OPTIONS: global
 } "Decrement loop counter"
 
-
+#########################################################################################
+# We must schedule the local routines to compute the integrals in global-loop-local 
+# (as opposed to local) mode, in order to make sure that the scheduling condition
+# "AFTER", which describes the ADMMass group, is respected. This group may depend 
+# on local routines (excision, emask) and must be run after all local routines. 
+# If we had no such a dependence, we could have scheduled the integral computations 
+# simply in local mode.
+##########################################################################################
 schedule GROUP ADMMass AT POSTSTEP AFTER ADMMass_InitLoopCounter WHILE ADMMass::ADMMass_LoopCounter
 {
   STORAGE:ADMMass_GFs[3]
   STORAGE:ADMMass_box
+  STORAGE:grid_spacing_product
 } "ADMMass loop"
 
 schedule ADMMass_Surface IN ADMMass
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]);
       }
-}
 
+}