- update to use new return codes from SpaceMask (-1=error instead of 0=error)

- split up a single huge function into several smaller/more readable
  subfunctions
- general cleanup of the logic for the noshrink option
- more documentation on same


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@1061 f88db872-0e4f-0410-b76b-b9085cfa78c5

1 files changed, 222 insertions, 128 deletions

diff --git a/src/driver/mask.cc b/src/driver/mask.cc
index e945832..ef963d4 100644
--- a/src/driver/mask.cc
+++ b/src/driver/mask.cc
@@ -4,6 +4,8 @@
 // <<<prototypes for functions local to this file>>>
 // setup_mask_dataptrs_and_bitfields - map gridfn/bitfield names to ptr/bitmasks
 // set_mask_gridfn - set mask gridfn(s) based on each horizon's shape
+/// set_mask_gridfn_to_outside_value - ... "outside" value
+/// set_mask_gridfn_to_inside_and_buffer_values - "inside"/"buffer" values
 //
 
 #include <stdio.h>
@@ -46,6 +48,26 @@ using jtutil::error_exit;
 //******************************************************************************
 
 //
+// prototypes for functions local to this file
+//
+
+namespace {
+void set_mask_gridfn_to_outside_value(const struct cactus_grid_info& cgi,
+				      const struct mask_info& mask_info,
+				      const struct verbose_info& verbose_info);
+void set_mask_gridfn_to_inside_and_buffer_values
+	(const struct cactus_grid_info& cgi,
+	 int use_min_i, int use_max_i,
+	 int use_min_j, int use_max_j,
+	 int use_min_k, int use_max_k,
+	 const struct mask_info& mask_info,
+	 const patch_system& ps,
+	 const struct verbose_info& verbose_info);
+	  }
+
+//******************************************************************************
+
+//
 // This function maps the character-string names of the mask gridfn(s)
 // and/or bitfield(s) into internal data pointers and/or bit masks/values:
 //
@@ -58,13 +80,10 @@ using jtutil::error_exit;
 //		 .buffer_value    --> .buffer_bitvalue
 //		 .outside_value   --> .outside_bitvalue
 //
-// FIXME: right now the SpaceMask() functions return 0 for error;
-//        this may change to -1 in the future
-//
 void setup_mask_dataptrs_and_bitfields(const cGH *GH,
 				       struct mask_info& mask_info)
 {
-const int SPACEMASK_ERROR = 0;
+const int SPACEMASK_ERROR = -1;
 
 struct mask_info::old_style_mask_info& osmi = mask_info.old_style_mask_info;
 if (mask_info.set_old_style_mask)
@@ -124,12 +143,6 @@ if (mask_info.set_new_style_mask)
 // This function sets the mask grid function specified by  mask_info
 // based on each horizon's shape.
 //
-// FIXME: The current implementation is quite inefficient: it does
-//        a full 2-D interpolation for each xyz grid point within
-//        each horizon's xyz bounding box.  It would be more efficient
-//        to batch the interpolations, and maybe also use r_min/r_max
-//        tests for early-out...
-//
 void set_mask_gridfn(int N_horizons,
 		     const struct AH_data* const AH_data_array[],
 		     const struct cactus_grid_info& cgi,
@@ -176,53 +189,13 @@ if (verbose_info.print_algorithm_debug)
 	}
 
 
-
 //
 // set the mask to the outside value everywhere in
 // (this processor's chunk of) the grid
 //
-if (verbose_info.print_algorithm_details)
-   then CCTK_VInfo(CCTK_THORNSTRING,
-		   "   setting mask grid function to \"outside\"");
-
-	for (int k = 0 ; k < cgi.local_gridfn_dims[Z_AXIS] ; ++k)
-	{
-	for (int j = 0 ; j < cgi.local_gridfn_dims[Y_AXIS] ; ++j)
-	{
-	for (int i = 0 ; i < cgi.local_gridfn_dims[X_AXIS] ; ++i)
-	{
-	const int posn = CCTK_GFINDEX3D(cgi.GH, i,j,k);
-	if (set_old_style_mask)
-	   then {
-		if (mask_info.mask_is_noshrink)
-		   then {
-			CCTK_REAL old_value = osmi.gridfn_dataptr[posn];
-			if ((old_value != osmi.inside_value) && (old_value != osmi.buffer_value))
-			   then osmi.gridfn_dataptr[posn] = osmi.outside_value;
-			}
-		   else osmi.gridfn_dataptr[posn] = osmi.buffer_value;
-		}
-	if (set_new_style_mask)
-	   then {
-		if (mask_info.mask_is_noshrink)
-		   then {
-			if (    !SpaceMask_CheckStateBits(nsmi.gridfn_dataptr, posn,
-							  nsmi.bitfield_bitmask,
-							  nsmi.inside_bitvalue)
-			     && !SpaceMask_CheckStateBits(nsmi.gridfn_dataptr, posn,
-							  nsmi.bitfield_bitmask,
-							  nsmi.buffer_bitvalue)    )
-			   then SpaceMask_SetStateBits(nsmi.gridfn_dataptr, posn,
-						       nsmi.bitfield_bitmask,
-						       nsmi.outside_bitvalue);
-			}
-		   else SpaceMask_SetStateBits(nsmi.gridfn_dataptr, posn,
-					       nsmi.bitfield_bitmask,
-					       nsmi.outside_bitvalue);
-		}
-	}
-	}
-	}
+set_mask_gridfn_to_outside_value(cgi,
+				 mask_info,
+				 verbose_info);
 
 
 //
@@ -244,6 +217,7 @@ if (verbose_info.print_algorithm_details)
 	const patch_system& ps                      = *AH_data.ps_ptr;
 	const struct BH_diagnostics& BH_diagnostics =  AH_data.BH_diagnostics;
 
+
 	// horizon bounding box, rounded "out" to the next grid point
 	const int AH_min_i = cgi.ijk_floor_of_global_xyz(X_AXIS, BH_diagnostics.min_x);
 	const int AH_max_i = cgi.ijk_ceil_of_global_xyz (X_AXIS, BH_diagnostics.max_x);
@@ -308,90 +282,210 @@ if (verbose_info.print_algorithm_details)
 			   double(cgi.global_xyz_of_ijk(Z_AXIS, use_max_k)));
 		}
 
-	//
-	// set the mask for each point in the intersection
-	//
-	long inside_count = 0;
-	long buffer_count = 0;
-		for (int k = use_min_k ; k <= use_max_k ; ++k)
-		{
-		for (int j = use_min_j ; j <= use_max_j ; ++j)
-		{
-		for (int i = use_min_i ; i <= use_max_i ; ++i)
-		{
-		const int posn = CCTK_GFINDEX3D(cgi.GH, i,j,k);
-
-		const fp global_x = cgi.global_xyz_of_ijk(X_AXIS, i);
-		const fp global_y = cgi.global_xyz_of_ijk(Y_AXIS, j);
-		const fp global_z = cgi.global_xyz_of_ijk(Z_AXIS, k);
-		const fp local_x = global_x - ps.origin_x();
-		const fp local_y = global_y - ps.origin_y();
-		const fp local_z = global_z - ps.origin_z();
-
-		const fp r = jtutil::hypot3(local_x, local_y, local_z);
-
-		// interpolate to find r_horizon in this direction
-		//
-		// FIXME:
-		// it would be more efficient here to compute the
-		// radia of a whole batch of points at once
-		const fp r_horizon = ps.radius_in_local_xyz_direction
-					   (gfns::gfn__h,
-					    local_x, local_y, local_z);
-
-		const fp r_inner = mask_info.radius_multiplier*r_horizon
-				   + mask_info.radius_offset*Cactus_dx;
-		const fp r_outer = r_inner + mask_info.buffer_thickness*Cactus_dx;
-
-		if	(r <= r_inner)
+	set_mask_gridfn_to_inside_and_buffer_values(cgi,
+						    use_min_i, use_max_i,
+						    use_min_j, use_max_j,
+						    use_min_k, use_max_k,
+						    mask_info,
+						    ps,
+						    verbose_info);
+	}
+}
+
+//******************************************************************************
+
+//
+// For each Cactus grid point in (this processor's chunk of) the Cactus
+// grid, this function sets the mask grid function specified by  mask_info
+// to the appropriate value for the "outside" region (possibly leaving
+// the mask unchanged at its previous value if the  mask_is_noshrink
+// option is set).
+//
+namespace {
+void set_mask_gridfn_to_outside_value(const struct cactus_grid_info& cgi,
+				      const struct mask_info& mask_info,
+				      const struct verbose_info& verbose_info)
+{
+const bool set_old_style_mask = mask_info.set_old_style_mask;
+const bool set_new_style_mask = mask_info.set_new_style_mask;
+const struct mask_info::old_style_mask_info& osmi = mask_info.old_style_mask_info;
+const struct mask_info::new_style_mask_info& nsmi = mask_info.new_style_mask_info;
+
+
+if (verbose_info.print_algorithm_details)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   setting mask grid function to \"outside\"");
+
+	for (int k = 0 ; k < cgi.local_gridfn_dims[Z_AXIS] ; ++k)
+	{
+	for (int j = 0 ; j < cgi.local_gridfn_dims[Y_AXIS] ; ++j)
+	{
+	for (int i = 0 ; i < cgi.local_gridfn_dims[X_AXIS] ; ++i)
+	{
+	const int posn = CCTK_GFINDEX3D(cgi.GH, i,j,k);
+	if (set_old_style_mask)
+	   then {
+		const CCTK_REAL old_value = osmi.gridfn_dataptr[posn];
+		if ( mask_info.mask_is_noshrink
+		     && (    (old_value == osmi.inside_value)
+			  || (old_value == osmi.buffer_value)    ) )
 		   then {
-			++inside_count;
-			if (set_old_style_mask)
-			   then osmi.gridfn_dataptr[posn] = osmi.inside_value;
-			if (set_new_style_mask)
-			   then SpaceMask_SetStateBits(nsmi.gridfn_dataptr, posn,
-						       nsmi.bitfield_bitmask,
-						       nsmi.inside_bitvalue);
+			// we're in "noshrink mode"
+			// and this point was previously "inside" or "buffer"
+			// ==> no-op here
 			}
-		else if (r <= r_outer)
+		   else osmi.gridfn_dataptr[posn] = osmi.buffer_value;
+		}
+	if (set_new_style_mask)
+	   then {
+		if ( mask_info.mask_is_noshrink
+		     && (    SpaceMask_CheckStateBits(nsmi.gridfn_dataptr, posn,
+						      nsmi.bitfield_bitmask,
+						      nsmi.inside_bitvalue)
+			  || SpaceMask_CheckStateBits(nsmi.gridfn_dataptr, posn,
+						      nsmi.bitfield_bitmask,
+						      nsmi.buffer_bitvalue)    ) )
 		   then {
-			++buffer_count;
-			if (set_old_style_mask)
+			// we're in "noshrink mode"
+			// and this point was previously "inside" or "buffer"
+			// ==> no-op here
+			}
+		   else SpaceMask_SetStateBits(nsmi.gridfn_dataptr, posn,
+					       nsmi.bitfield_bitmask,
+					       nsmi.outside_bitvalue);
+		}
+	}
+	}
+	}
+}
+	  }
+
+//******************************************************************************
+
+//
+// For each Cactus grid point in a specified region of the Cactus grid,
+// and for a specified (single) horizon, this function either sets the
+// mask grid function specified by  mask_info to the appropriate value
+// if the point is in the "inside" or "buffer" region for a specified
+// horizon (possibly leaving the mask unchanged at its previous value
+// if the  mask_is_noshrink  option is set), or leaves the mask grid
+// function unchanged if the point is actually outside this horizon.
+//
+// Arguments:
+// ps = The patch system for this horizon.
+// use_{min,max}_{i,j,k} = Specifies the region of the Cactus grid.
+//
+// FIXME:
+// The current implementation is quite inefficient: it does a full 2-D
+// interpolation (to find the horizon radius) for each xyz grid point
+// in the specified region.  It would be more efficient to batch the
+// interpolations, and maybe also use r_min/r_max tests for early-out.
+//
+namespace {
+void set_mask_gridfn_to_inside_and_buffer_values
+	(const struct cactus_grid_info& cgi,
+	 int use_min_i, int use_max_i,
+	 int use_min_j, int use_max_j,
+	 int use_min_k, int use_max_k,
+	 const struct mask_info& mask_info,
+	 const patch_system& ps,
+	 const struct verbose_info& verbose_info)
+{
+const bool set_old_style_mask = mask_info.set_old_style_mask;
+const bool set_new_style_mask = mask_info.set_new_style_mask;
+const struct mask_info::old_style_mask_info& osmi = mask_info.old_style_mask_info;
+const struct mask_info::new_style_mask_info& nsmi = mask_info.new_style_mask_info;
+const fp Cactus_dx = cgi.mean_coord_delta;
+
+long inside_count = 0;
+long buffer_count = 0;
+
+	for (int k = use_min_k ; k <= use_max_k ; ++k)
+	{
+	for (int j = use_min_j ; j <= use_max_j ; ++j)
+	{
+	for (int i = use_min_i ; i <= use_max_i ; ++i)
+	{
+	const int posn = CCTK_GFINDEX3D(cgi.GH, i,j,k);
+
+	const fp global_x = cgi.global_xyz_of_ijk(X_AXIS, i);
+	const fp global_y = cgi.global_xyz_of_ijk(Y_AXIS, j);
+	const fp global_z = cgi.global_xyz_of_ijk(Z_AXIS, k);
+	const fp local_x = global_x - ps.origin_x();
+	const fp local_y = global_y - ps.origin_y();
+	const fp local_z = global_z - ps.origin_z();
+
+	const fp r = jtutil::hypot3(local_x, local_y, local_z);
+
+	// interpolate to find r_horizon in this direction
+	//
+	// FIXME:
+	// it would be more efficient here to compute the
+	// radia of a whole batch of points at once
+	const fp r_horizon = ps.radius_in_local_xyz_direction
+				   (gfns::gfn__h,
+				    local_x, local_y, local_z);
+
+	const fp r_inner = mask_info.radius_multiplier*r_horizon
+			   + mask_info.radius_offset*Cactus_dx;
+	const fp r_outer = r_inner + mask_info.buffer_thickness*Cactus_dx;
+
+	if	(r <= r_inner)
+	   then {
+		// set the mask to the "inside" value
+		++inside_count;
+		if (set_old_style_mask)
+		   then osmi.gridfn_dataptr[posn] = osmi.inside_value;
+		if (set_new_style_mask)
+		   then SpaceMask_SetStateBits(nsmi.gridfn_dataptr, posn,
+					       nsmi.bitfield_bitmask,
+					       nsmi.inside_bitvalue);
+		}
+	else if (r <= r_outer)
+	   then {
+		// set the mask to the "buffer" value
+		// ... except that if the  mask_is_noshrink  option is set
+		//     and the mask is already set to the "inside" value
+		//     then we need to leave it unchanged
+		//     (at the "inside" value)
+		++buffer_count;
+		if (set_old_style_mask)
+		   then {
+			if ( mask_info.mask_is_noshrink
+			     && (osmi.gridfn_dataptr[posn] == osmi.inside_value) )
 			   then {
-				if (mask_info.mask_is_noshrink)
-				   then {
-					if (osmi.gridfn_dataptr[posn] != osmi.inside_value)
-					   then osmi.gridfn_dataptr[posn] = osmi.buffer_value;
-					}
-				   else osmi.gridfn_dataptr[posn] = osmi.buffer_value;
+				// we're in "noshrink mode"
+				// and this point was previously "inside"
+				// ==> no-op here
 				}
-			if (set_new_style_mask)
+			   else osmi.gridfn_dataptr[posn] = osmi.buffer_value;
+			}
+		if (set_new_style_mask)
+		   then {
+			if ( mask_info.mask_is_noshrink
+			     && SpaceMask_CheckStateBits(nsmi.gridfn_dataptr, posn,
+							 nsmi.bitfield_bitmask,
+							 nsmi.inside_bitvalue) )
 			   then {
-				if (mask_info.mask_is_noshrink)
-				   then {
-					if (! SpaceMask_CheckStateBits(nsmi.gridfn_dataptr, posn,
-								       nsmi.bitfield_bitmask,
-								       nsmi.inside_bitvalue))
-					   then SpaceMask_SetStateBits(nsmi.gridfn_dataptr, posn,
-								       nsmi.bitfield_bitmask,
-								       nsmi.buffer_bitvalue);
-					}
-				   else SpaceMask_SetStateBits(nsmi.gridfn_dataptr, posn,
-							       nsmi.bitfield_bitmask,
-							       nsmi.buffer_bitvalue);
+				// we're in "noshrink mode"
+				// and this point was previously "inside"
+				// ==> no-op here
 				}
-			}
-		else	{
-			// it's "outside", and we've already done that case above
-			// ==> no-op here
+			   else SpaceMask_SetStateBits(nsmi.gridfn_dataptr, posn,
+						       nsmi.bitfield_bitmask,
+						       nsmi.buffer_bitvalue);
 			}
 		}
+	else	{
+		// it's "outside" ==> no-op here
 		}
-		}
-
-	if (verbose_info.print_algorithm_details)
-	   then CCTK_VInfo(CCTK_THORNSTRING,
-			   "      %ld \"inside\" points, %ld \"buffer\" points on this processor",
-			   inside_count, buffer_count);
 	}
+	}
+	}
+
+if (verbose_info.print_algorithm_details)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "      %ld \"inside\" points, %ld \"buffer\" points on this processor",
+		   inside_count, buffer_count);
 }
+	  }