new file I forgot to add before -- needed for the

(experimental, not-yet-tested-well, not-yet-documented-at-all)
function-aliasing interface for other thorns to query the AH shape


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

1 files changed, 145 insertions, 0 deletions

diff --git a/src/driver/horizon_radius.cc b/src/driver/horizon_radius.cc
new file mode 100644
index 0000000..a5982b6
--- /dev/null
+++ b/src/driver/horizon_radius.cc
@@ -0,0 +1,145 @@
+// horizon_radius.cc -- provide horizon radius & other info for other thorns
+// $Header$
+//
+// <<<access to persistent data>>>
+// AHFinderDirect_local_coordinate_origin - provide our local coordinate origin
+// AHFinderDirect_radius_in_direction - provide r(angle) function
+//
+
+#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 "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"
+
+// all the code in this file is inside this namespace
+namespace AHFinderDirect
+	  {
+
+//******************************************************************************
+
+//
+// ***** access to persistent data *****
+//
+extern struct state state;
+
+//******************************************************************************
+
+//
+// This function is called (via the magic of function aliasing) by
+// other thorns to find out our local coordinate origin for a given AH.
+//
+extern "C"
+  void AHFinderDirect_local_coordinate_origin
+     (CCTK_INT horizon_number,
+      CCTK_REAL* origin_x_ptr, CCTK_REAL* origin_y_ptr, CCTK_REAL* origin_z_ptr)
+{
+const struct verbose_info& verbose_info = state.verbose_info;
+
+if (!  ((horizon_number >= 1) && (horizon_number < state.N_horizons))  )
+   then CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"AHFinderDirect_local_coordinate_origin():\n"
+"        horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
+		   ,
+		   int(horizon_number), state.N_horizons);	/*NOTREACHED*/
+
+assert(state.AH_data_array[horizon_number] != NULL);
+const struct AH_data& AH_data = *state.AH_data_array[horizon_number];
+
+assert(AH_data.ps_ptr != NULL);
+const patch_system& ps = *AH_data.ps_ptr;
+
+*origin_x_ptr = ps.origin_x();
+*origin_y_ptr = ps.origin_y();
+*origin_z_ptr = ps.origin_z();
+}
+
+//******************************************************************************
+
+//
+// This function is called (via the magic of function aliasing) by
+// other thorns to find out a given AH's radius in the direction from
+// its local coordinate origin to a given (x,y,z) coordinate or coordinates.
+//
+// FIXME:
+// The current implementation is quite inefficient: it does a full 2-D
+// interpolation (to find the horizon radius) for each of the xyz points.
+// It would be more efficient to batch the interpolations.
+//
+// Arguments:
+// horizon_number = must be in the range 1 to N_horizons
+// N_points = should be >= 0
+// x[], y[], z[] = these give the (x,y,z) coordinates
+// radius[] = this is set to the horizon radius values (Euclidean distance
+//	      from the local coordinate origin), or to all -1.0 if we didn't
+//	      find this horizon the last time we looked for it
+//
+extern "C"
+  void AHFinderDirect_radius_in_direction
+	(CCTK_INT horizon_number,
+	 CCTK_INT N_points,
+	 const CCTK_REAL* const x, const CCTK_REAL* const y, const CCTK_REAL* const z,
+	 CCTK_REAL* const radius)
+{
+const struct verbose_info& verbose_info = state.verbose_info;
+
+if (! ((horizon_number >= 1) && (horizon_number < state.N_horizons)) )
+   then CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"AHFinderDirect_distance_outside_thorn():\n"
+"        horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
+		   ,
+		   int(horizon_number), state.N_horizons);	/*NOTREACHED*/
+
+assert(state.AH_data_array[horizon_number] != NULL);
+const struct AH_data& AH_data = *state.AH_data_array[horizon_number];
+
+	for (int point = 0 ; point < N_points ; ++point)
+	{
+	if (AH_data.found_flag)
+	   then {
+		assert(AH_data.ps_ptr != NULL);
+		const patch_system& ps = *AH_data.ps_ptr;
+
+		const fp local_x = x[point] - ps.origin_x();
+		const fp local_y = y[point] - ps.origin_y();
+		const fp local_z = z[point] - ps.origin_z();
+		radius[point] = ps.radius_in_local_xyz_direction
+				      (gfns::gfn__h,
+				       local_x, local_y, local_z);
+		}
+	   else radius[point] = -1.0;
+	}
+}
+
+//******************************************************************************
+
+	  }	// namespace AHFinderDirect