bring cvs back into line with files from my laptop

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

1 files changed, 337 insertions, 0 deletions

diff --git a/src/driver/mask.cc b/src/driver/mask.cc
new file mode 100644
index 0000000..8e66f1d
--- /dev/null
+++ b/src/driver/mask.cc
@@ -0,0 +1,337 @@
+// mask.cc -- set a mask gridfn based on each horizon shape
+// $Header$
+//
+// <<<prototypes for functions local to this file>>>
+// set_mask_gridfn - set mask gridfn based on each horizon's shape
+//
+
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+
+#include "util_Table.h"
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "SpaceMask.h"			// from thorn SpaceMask
+
+#include "config.h"
+#include "stdc.h"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
+using jtutil::error_exit;
+
+#include "../patch/coords.hh"
+#include "../patch/grid.hh"
+#include "../patch/fd_grid.hh"
+#include "../patch/patch.hh"
+#include "../patch/patch_edge.hh"
+#include "../patch/patch_interp.hh"
+#include "../patch/ghost_zone.hh"
+#include "../patch/patch_system.hh"
+
+#include "../elliptic/Jacobian.hh"
+
+#include "../gr/gfns.hh"
+#include "../gr/gr.hh"
+
+#include "horizon_sequence.hh"
+#include "BH_diagnostics.hh"
+#include "driver.hh"
+
+//******************************************************************************
+
+//
+// This function sets the mask grid function specified by  mask_info
+// based on each horizon's shape.
+//
+void set_mask_gridfn(int N_horizons,
+		     const struct AH_data* const AH_data_array[],
+		     const struct cactus_grid_info& cgi,
+		     const struct mask_info& mask_info,
+		     const struct verbose_info& verbose_info)
+{
+const struct mask_info::old_style_mask_info& old_style_mask_info
+	= mask_info.old_style_mask_info;
+const struct mask_info::new_style_mask_info& new_style_mask_info
+	= mask_info.new_style_mask_info;
+
+if (verbose_info.print_algorithm_highlights)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "setting mask grid function %s",
+		   mask_info.mask_gridfn_name);
+
+if (verbose_info.print_algorithm_debug)
+   then {
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "      grid on this processor has x=[%g,%g]",
+		   double(cgi.global_xyz_of_ijk(X_AXIS, 0)),
+		   double(cgi.global_xyz_of_ijk(X_AXIS, cgi.local_gridfn_dims[X_AXIS]-1)));
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "                                 y=[%g,%g]",
+		   double(cgi.global_xyz_of_ijk(Y_AXIS, 0)),
+		   double(cgi.global_xyz_of_ijk(Y_AXIS, cgi.local_gridfn_dims[Y_AXIS]-1)));
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "                                 z=[%g,%g]",
+		   double(cgi.global_xyz_of_ijk(Z_AXIS, 0)),
+		   double(cgi.global_xyz_of_ijk(Z_AXIS, cgi.local_gridfn_dims[Z_AXIS]-1)));
+	}
+
+
+//
+// set up the gridfn data pointers
+//
+CCTK_REAL* old_style_mask_gridfn_data_ptr = NULL;
+CCTK_INT * new_style_mask_gridfn_data_ptr = NULL;
+switch	(mask_info.mask_type)
+	{
+case mask_type__old_style:
+	old_style_mask_gridfn_data_ptr
+		= Cactus_gridfn_data_ptr<CCTK_REAL>
+					(cgi.GH,
+					 mask_info.mask_gridfn_varindex,
+					 mask_info.mask_gridfn_name);
+	break;
+case mask_type__new_style:
+	new_style_mask_gridfn_data_ptr
+		= Cactus_gridfn_data_ptr<CCTK_INT>
+					(cgi.GH,
+					 mask_info.mask_gridfn_varindex,
+					 mask_info.mask_gridfn_name);
+	break;
+default:
+	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   set_mask_gridfn(): unknown mask_info.mask_type=(int)%d\n"
+"                      when precomputing gridfn data pointers!\n"
+"                      (this should never happen!)"
+		   ,
+		   int(mask_info.mask_type));		/*NOTREACHED*/
+	}
+
+
+//
+// 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);
+	switch	(mask_info.mask_type)
+		{
+	case mask_type__old_style:
+		old_style_mask_gridfn_data_ptr[posn] = old_style_mask_info
+						       .outside_value;
+		break;
+	case mask_type__new_style:
+		SpaceMask_SetStateBits(new_style_mask_gridfn_data_ptr, posn,
+				       new_style_mask_info.bitfield_mask,
+				       new_style_mask_info.inside_value);
+		break;
+	CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   set_mask_gridfn(): unknown mask_info.mask_type=(int)%d\n"
+"                      when setting mask to outside value!\n"
+"                      (this should never happen!)"
+		   ,
+		   int(mask_info.mask_type));			/*NOTREACHED*/
+		}
+	}
+	}
+	}
+
+
+//
+// loop over each horizon's xyz bounding box
+// (intersected with (this processor's chunk of) the grid)
+// to set the mask accurately
+//
+	for (int hn = 1 ; hn <= N_horizons ; ++hn)
+	{
+	const struct AH_data& AH_data = *AH_data_array[hn];
+	if (! AH_data.found_flag)
+	   then continue;			// *** LOOP CONTROL ***
+
+	if (verbose_info.print_algorithm_details)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   setting mask grid function to \"inside\" for horizon %d",
+			   hn);
+
+	const patch_system& ps                      = *AH_data.ps_ptr;
+	const struct BH_diagnostics& BH_diagnostics =  AH_data.BH_diagnostics;
+
+	const fp radius_multiplier = mask_info.radius_multiplier[hn];
+	const fp radius_offset     = mask_info.radius_offset[hn]
+				     * cgi.mean_coord_delta;
+
+	// 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);
+	const int AH_min_j = cgi.ijk_floor_of_global_xyz(Y_AXIS, BH_diagnostics.min_y);
+	const int AH_max_j = cgi.ijk_ceil_of_global_xyz (Y_AXIS, BH_diagnostics.max_y);
+	const int AH_min_k = cgi.ijk_floor_of_global_xyz(Z_AXIS, BH_diagnostics.min_z);
+	const int AH_max_k = cgi.ijk_ceil_of_global_xyz (Z_AXIS, BH_diagnostics.max_z);
+	if (verbose_info.print_algorithm_debug)
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "      horizon bounding box is x=[%g,%g]",
+			   double(BH_diagnostics.min_x),
+			   double(BH_diagnostics.max_x));
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "                              y=[%g,%g]",
+			   double(BH_diagnostics.min_y),
+			   double(BH_diagnostics.max_y));
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "                              z=[%g,%g]",
+			   double(BH_diagnostics.min_z),
+			   double(BH_diagnostics.max_z));
+		}
+
+	//
+	// compute intersection of horizon bounding box with
+	// (this processor's chunk of) the grid
+	//
+	// as can be seen from the following picture, the intersection
+	// has the max/min of the min/max coordinate as its min/max:
+	//
+	//	           +--------------------+
+	//	           |                    |
+	//	+----------+--------------+     |
+	//	|          |XXXXXXXXXXXXXX|     |
+	//	|          |XXXXXXXXXXXXXX|     |
+	//	|          |XXXXXXXXXXXXXX|     |
+	//	|          +--------------+-----+
+	//	|                         |
+	//	|                         |
+	//	+-------------------------+
+	//
+	const int use_min_i = jtutil::max(AH_min_i, 0);
+	const int use_max_i = jtutil::min(AH_max_i, cgi.local_gridfn_dims[X_AXIS]-1);
+	const int use_min_j = jtutil::max(AH_min_j, 0);
+	const int use_max_j = jtutil::min(AH_max_j, cgi.local_gridfn_dims[Y_AXIS]-1);
+	const int use_min_k = jtutil::max(AH_min_k, 0);
+	const int use_max_k = jtutil::min(AH_max_k, cgi.local_gridfn_dims[Z_AXIS]-1);
+
+	if (verbose_info.print_algorithm_debug)
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "      use bounding box is x=[%g,%g]",
+			   double(cgi.global_xyz_of_ijk(X_AXIS, use_min_i)),
+			   double(cgi.global_xyz_of_ijk(X_AXIS, use_max_i)));
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "                          y=[%g,%g]",
+			   double(cgi.global_xyz_of_ijk(Y_AXIS, use_min_j)),
+			   double(cgi.global_xyz_of_ijk(Y_AXIS, use_max_j)));
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "                          z=[%g,%g]",
+			   double(cgi.global_xyz_of_ijk(Z_AXIS, use_min_k)),
+			   double(cgi.global_xyz_of_ijk(Z_AXIS, use_max_k)));
+		}
+
+	//
+	// set the mask for each point in the intersection
+	//
+	long inside_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);
+
+		CCTK_REAL old_style_mask_value;
+		CCTK_INT  new_style_mask_value;
+		if	(r > BH_diagnostics.max_radius)
+		   then {
+			// definitely outside
+			old_style_mask_value = old_style_mask_info.outside_value;
+			new_style_mask_value = new_style_mask_info.outside_value;
+			}
+		else if (r < BH_diagnostics.min_radius)
+		   then {
+			// definitely inside
+			++inside_count;
+			old_style_mask_value = old_style_mask_info.inside_value;
+			new_style_mask_value = new_style_mask_info.inside_value;
+			}
+		else	{
+			// interpolate to find out inside/outside status
+			//
+			// 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_mask = radius_multiplier*r_horizon
+					  + radius_offset;
+			if (r < r_mask)
+			   then {
+				++inside_count;
+				old_style_mask_value
+					= old_style_mask_info.inside_value;
+				new_style_mask_value
+					= new_style_mask_info.inside_value;
+				}
+			   else {
+				old_style_mask_value
+					= old_style_mask_info.outside_value;
+				new_style_mask_value
+					= new_style_mask_info.outside_value;
+				}
+			}
+
+		// set the mask
+		switch	(mask_info.mask_type)
+			{
+		case mask_type__old_style:
+			old_style_mask_gridfn_data_ptr[posn]
+				= old_style_mask_value;
+			break;
+		case mask_type__new_style:
+			SpaceMask_SetStateBits(new_style_mask_gridfn_data_ptr,
+					       posn,
+					       new_style_mask_info.bitfield_mask,
+					       new_style_mask_value);
+			break;
+		default:
+			CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__,
+				   CCTK_THORNSTRING,
+"\n"
+"   set_mask_gridfn(): unknown mask_info.mask_type=(int)%d!\n"
+"                      when setting value for horizon!\n"
+"                      (this should never happen!)"
+				   ,
+				   int(mask_info.mask_type));	/*NOTREACHED*/
+			}
+		}
+		}
+		}
+
+	if (verbose_info.print_algorithm_details)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "   ==> %ld points set to \"inside\"",
+			   inside_count);
+	}
+}