fix syntax errors in the last set of checkins

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

2 files changed, 171 insertions, 39 deletions

diff --git a/src/gr/driver.cc b/src/gr/driver.cc
index 58eeaf6..248797c 100644
--- a/src/gr/driver.cc
+++ b/src/gr/driver.cc
@@ -1,7 +1,9 @@
 // driver.cc -- top level driver for finding apparent horizons
 // $Id$
 //
+// <<<prototypes for functions local to this file>>>
 // AHFinderDirect_driver - top-level driver
+/// setup_ellipsoid_initial_guess - set up ellipsoid in h gridfn
 //
 
 #include <stdio.h>
@@ -36,7 +38,18 @@ using jtutil::error_exit;
 #include "AHFinderDirect.hh"
 
 //******************************************************************************
-//******************************************************************************
+
+//
+// ***** prototypes for functions local to this file *****
+//
+
+namespace {
+void setup_ellipsoid_initial_guess
+	(patch_system& ps,
+	 fp global_center_x, fp global_center_y, fp global_center_z,
+	 fp radius_x, fp radius_y, fp radius_z);
+	  };
+
 //******************************************************************************
 
 //
@@ -50,6 +63,25 @@ DECLARE_CCTK_PARAMETERS
 
 CCTK_VInfo(CCTK_THORNSTRING, "initializing AHFinderDirect data structures");
 
+
+//
+// set up the geometry interpolator
+//
+struct geometry_interpolator_info gii;
+CCTK_VInfo(CCTK_THORNSTRING, "   setting up geometry interpolator");
+gii.operator_handle = CCTK_InterpHandle(geometry_interpolator_name);
+if (gii.operator_handle < 0)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "couldn't find interpolator \"%s\"!",
+		   geometry_interpolator_name);		/*NOTREACHED*/
+
+gii.param_table_handle = Util_TableCreateFromString(geometry_interpolator_pars);
+if (gii.param_table_handle < 0)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "bad geometry-interpolator parameter(s) \"%s\"!",
+		   geometry_interpolator_pars);		/*NOTREACHED*/
+
+
 //
 // set up the interpatch interpolator
 //
@@ -59,13 +91,14 @@ if (interp_handle < 0)
    then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "couldn't find interpolator \"%s\"!",
 		   interpatch_interpolator_name);		/*NOTREACHED*/
-const int interp_par_table_handle
+const int interp_param_table_handle
 	= Util_TableCreateFromString(interpatch_interpolator_pars);
-if (interp_par_table_handle < 0)
+if (interp_param_table_handle < 0)
    then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "bad interpatch-interpolator parameter(s) \"%s\"!",
 		   interpatch_interpolator_pars);		/*NOTREACHED*/
 
+
 //
 // set up the Cactus grid info
 //
@@ -118,22 +151,6 @@ cgi.K_dd_33_data = static_cast<const fp *>(
 			CCTK_VarDataPtr(cctkGH, 0, "einstein::kzz")
 					  );
 
-//
-// set up the geometry interpolator
-//
-struct geometry_interpolator_info gii;
-CCTK_VInfo(CCTK_THORNSTRING, "   setting up geometry interpolator");
-gii.operator_handle = CCTK_InterpHandle(geometry_interpolator_name);
-if (interp_handle < 0)
-   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "couldn't find interpolator \"%s\"!",
-		   geometry_interpolator_name);		/*NOTREACHED*/
-gii.param_table_handle
-	= Util_TableCreateFromString(geometry_interpolator_pars);
-if (interp_par_table_handle < 0)
-   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		   "bad geometry-interpolator parameter(s) \"%s\"!",
-		   geometry_interpolator_pars);		/*NOTREACHED*/
 
 //
 // create the patch system and initialize the xyz derivative coefficients
@@ -144,14 +161,39 @@ patch_system ps(origin_x, origin_y, origin_z,
 		N_ghost_points, N_overlap_points, delta_drho_dsigma,
 		nominal_gfns::min_gfn, nominal_gfns::max_gfn,
 		ghosted_gfns::min_gfn, ghosted_gfns::max_gfn,
-		interp_handle, interp_par_table_handle);
+		interp_handle, interp_param_table_handle);
+
+
+//
+// set up the initial guess for the apparent horizon shape
+//
+if	(STRING_EQUAL(initial_guess_method, "read from file"))
+   then {
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   reading initial guess from \"%s\"",
+		   initial_guess__read_from_file__file_name);
+	ps.read_ghosted_gridfn(ghosted_gfns::gfn__h,
+			       initial_guess__read_from_file__file_name,
+			       false);		// no ghost zones
+	}
+else if (STRING_EQUAL(initial_guess_method, "ellipsoid"))
+   then setup_ellipsoid_initial_guess(ps,
+				      initial_guess__ellipsoid__center_global_x,
+				      initial_guess__ellipsoid__center_global_y,
+				      initial_guess__ellipsoid__center_global_z,
+				      initial_guess__ellipsoid__radius_x,
+				      initial_guess__ellipsoid__radius_y,
+				      initial_guess__ellipsoid__radius_z);
+else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "unknown initial_guess_method=\"%s\"!",
+		   initial_guess_method);			/*NOTREACHED*/
+
 
 //
 // find the apparent horizon
 //
 if      (STRING_EQUAL(method, "horizon"))
    then {
-	ps.read_ghosted_gridfn(ghosted_gfns::gfn__h, "h.dat", false);
 	horizon_function(ps, cgi, gii);
 	ps.print_gridfn(nominal_gfns::gfn__H, "H.dat");
 	}
@@ -159,3 +201,93 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "unknown method=\"%s\"!",
 		   method);					/*NOTREACHED*/
 }
+
+//******************************************************************************
+
+//
+// This function sets up an ellipsoid in the gridfn h, using the
+// formulas in "ellipsoid.maple" and the Maple-generated C code in
+// "ellipsoid.c":
+//
+// ellipsoid has center (A,B,C), radius (a,b,c)
+// angular coordinate system has center (U,V,W)
+//
+// direction cosines wrt angular coordinate center are (alpha,beta,gamma)
+// i.e. a point has coordinates (U+alpha*r, V+beta*r, W+gamma*r)
+//
+// then the equation of the ellipsoid is
+//	(U+alpha*r - A)^2     (V+beta*r - B)^2     (W+gamma*r - C)^2
+//	-----------------  +  ----------------  +  -----------------  =  1
+//	        a^2                  b^2                   c^2
+//
+// to solve this, we introduce intermediate variables
+//	AU = A - U
+//	BV = B - V
+//	CW = C - W
+//
+namespace {
+void setup_ellipsoid_initial_guess
+	(patch_system& ps,
+	 fp global_center_x, fp global_center_y, fp global_center_z,
+	 fp radius_x, fp radius_y, fp radius_z)
+{
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "   h = ellipsoid: global_center=(%g,%g,%g)",
+	   global_center_x, global_center_y, global_center_z);
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "                  radius=(%g,%g,%g)",
+	   radius_x, radius_y, radius_z);
+
+	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
+	{
+	patch& p = ps.ith_patch(pn);
+
+		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
+		{
+		for (int isigma = p.min_isigma() ;
+		     isigma <= p.max_isigma() ;
+		     ++isigma)
+		{
+		const fp rho = p.rho_of_irho(irho);
+		const fp sigma = p.sigma_of_isigma(isigma);
+		fp alpha, beta, gamma;
+		p.alpha_beta_gamma_of_rho_sigma(rho,sigma, alpha,beta,gamma);
+
+		// set up variables used by Maple-generated code
+		const fp AU = global_center_x - ps.origin_x();
+		const fp BV = global_center_y - ps.origin_y();
+		const fp CW = global_center_z - ps.origin_z();
+		const fp a = radius_x;
+		const fp b = radius_y;
+		const fp c = radius_z;
+
+		// compute the solutions r_plus and r_minus
+		fp r_plus, r_minus;
+		#include "ellipsoid.c"
+
+		// exactly one of the solutions (=def r) should be positive
+		fp r;
+		if      ((r_plus > 0.0) && (r_minus < 0.0))
+		   then r = r_plus;
+		else if ((r_plus < 0.0) && (r_minus > 0.0))
+		   then r = r_minus;
+		else    CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+				   "\n"
+"   expected exactly one r>0 solution, got 0 or 2!\n"
+"   %s patch (irho,isigma)=(%d,%d) ==> (rho,sigma)=(%g,%g)\n"
+"   (alpha,beta,gamma)=(%g,%g,%g)\n"
+"   ==> r_plus=%g r_minus=%g\n"
+				   ,
+				   p.name(), irho, isigma,
+				   double(rho), double(sigma),
+				   double(alpha), double(beta), double(gamma),
+				   double(r_plus), double(r_minus));
+		   						/*NOTREACHED*/
+
+		// r = horizon radius at this grid point
+		p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) = r;
+		}
+		}
+	}
+}
+	  }
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index ed1ddae..0cb1f7e 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -143,43 +143,43 @@ CCTK_VInfo(CCTK_THORNSTRING, "      xyz positions and derivative coefficients");
 
 		// 1st derivative coefficient gridfns X_ud
 		p.gridfn(nominal_gfns::gfn__X_ud_11, irho,isigma)
-			= p.partial_rho_wrt_x(x,y,z);
+			= p.partial_rho_wrt_x(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_ud_12, irho,isigma)
-			= p.partial_rho_wrt_y(x,y,z);
+			= p.partial_rho_wrt_y(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_ud_13, irho,isigma)
-			= p.partial_rho_wrt_z(x,y,z);
+			= p.partial_rho_wrt_z(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_ud_21, irho,isigma)
-			= p.partial_sigma_wrt_x(x,y,z);
+			= p.partial_sigma_wrt_x(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_ud_22, irho,isigma)
-			= p.partial_sigma_wrt_y(x,y,z);
+			= p.partial_sigma_wrt_y(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_ud_23, irho,isigma)
-			= p.partial_sigma_wrt_z(x,y,z);
+			= p.partial_sigma_wrt_z(local_x,local_y,local_z);
 
 		// 2nd derivative coefficient gridfns X_udd
 		p.gridfn(nominal_gfns::gfn__X_udd_111, irho,isigma)
-			= p.partial2_rho_wrt_xx(x,y,z);
+			= p.partial2_rho_wrt_xx(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_udd_112, irho,isigma)
-			= p.partial2_rho_wrt_xy(x,y,z);
+			= p.partial2_rho_wrt_xy(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_udd_113, irho,isigma)
-			= p.partial2_rho_wrt_xz(x,y,z);
+			= p.partial2_rho_wrt_xz(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_udd_122, irho,isigma)
-			= p.partial2_rho_wrt_yy(x,y,z);
+			= p.partial2_rho_wrt_yy(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_udd_123, irho,isigma)
-			= p.partial2_rho_wrt_yz(x,y,z);
+			= p.partial2_rho_wrt_yz(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_udd_133, irho,isigma)
-			= p.partial2_rho_wrt_zz(x,y,z);
+			= p.partial2_rho_wrt_zz(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_udd_211, irho,isigma)
-			= p.partial2_sigma_wrt_xx(x,y,z);
+			= p.partial2_sigma_wrt_xx(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_udd_212, irho,isigma)
-			= p.partial2_sigma_wrt_xy(x,y,z);
+			= p.partial2_sigma_wrt_xy(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_udd_213, irho,isigma)
-			= p.partial2_sigma_wrt_xz(x,y,z);
+			= p.partial2_sigma_wrt_xz(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_udd_222, irho,isigma)
-			= p.partial2_sigma_wrt_yy(x,y,z);
+			= p.partial2_sigma_wrt_yy(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_udd_223, irho,isigma)
-			= p.partial2_sigma_wrt_yz(x,y,z);
+			= p.partial2_sigma_wrt_yz(local_x,local_y,local_z);
 		p.gridfn(nominal_gfns::gfn__X_udd_233, irho,isigma)
-			= p.partial2_sigma_wrt_zz(x,y,z);
+			= p.partial2_sigma_wrt_zz(local_x,local_y,local_z);
 		}
 		}
 	}