change from $Id:$ to $Header:$

fix a bug in local_coords::partial_xyz_wrt_mu_nu()
redo Jacobian computation for surface integrals


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

23 files changed, 171 insertions, 66 deletions

diff --git a/src/patch/coord_derivs.maple b/src/patch/coord_derivs.maple
index 154eb20..afddb02 100644
--- a/src/patch/coord_derivs.maple
+++ b/src/patch/coord_derivs.maple
@@ -1,5 +1,5 @@
 # coord_derivs.maple -- compute derivatives of (mu,nu,phi) wrt (x,y,z)
-# $Id$
+# $Header$
 
 mu  := arctan(y/z);
 nu  := arctan(x/z);
diff --git a/src/patch/coords.cc b/src/patch/coords.cc
index 67c41f9..fefd2f1 100644
--- a/src/patch/coords.cc
+++ b/src/patch/coords.cc
@@ -1,5 +1,5 @@
 // coords.cc - coordinate conversions etc
-// $Id$
+// $Header$
 
 //
 // local_coords::fuzzy_EQ_{ang,dang}
@@ -294,7 +294,7 @@ partial_z_wrt_mu = -0.5 * (r2/(z*t2)) * partial_t_wrt_mu;
 partial_z_wrt_nu = -0.5 * (r2/(z*t2)) * partial_t_wrt_nu;
 
 partial_x_wrt_mu = tan_nu*partial_z_wrt_mu;
-partial_x_wrt_nu = tan_nu*partial_z_wrt_nu + z*(1.0+tan2_mu);
+partial_x_wrt_nu = tan_nu*partial_z_wrt_nu + z*(1.0+tan2_nu);
 partial_y_wrt_mu = tan_mu*partial_z_wrt_mu + z*(1.0+tan2_mu);
 partial_y_wrt_nu = tan_mu*partial_z_wrt_nu;
 }
diff --git a/src/patch/coords.hh b/src/patch/coords.hh
index bad88ea..787c203 100644
--- a/src/patch/coords.hh
+++ b/src/patch/coords.hh
@@ -1,5 +1,5 @@
 // coords.hh -- coordinate systems and conversions
-// $Id$
+// $Header$
 //
 // ***** coordinate systems *****
 // ***** conversions between different local coordinate systems *****
diff --git a/src/patch/fd_grid.cc b/src/patch/fd_grid.cc
index 5719316..d2b4638 100644
--- a/src/patch/fd_grid.cc
+++ b/src/patch/fd_grid.cc
@@ -1,5 +1,5 @@
 // fd_grid.cc -- grid with finite differencing operations
-// $Id$
+// $Header$
 //
 // fd_grid::dx_coeff
 // fd_grid::dxx_coeff
diff --git a/src/patch/fd_grid.hh b/src/patch/fd_grid.hh
index d0f511e..78ae6b6 100644
--- a/src/patch/fd_grid.hh
+++ b/src/patch/fd_grid.hh
@@ -1,5 +1,5 @@
 // fd_grid.hh -- grid with finite differencing operations
-// $Id$
+// $Header$
 //
 // *** Implementation Notes -- Overview ***
 // *** Implementation Notes -- Techniques using C++ Templates ***
diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc
index 3ace6e3..747cdb0 100644
--- a/src/patch/ghost_zone.cc
+++ b/src/patch/ghost_zone.cc
@@ -1,5 +1,5 @@
 // ghost_zone.cc -- fill in gridfn data in patch ghost zones
-// $Id$
+// $Header$
 
 //
 // ghost_zone::cast_to_symmetry_ghost_zone
diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh
index cb9e328..c548145 100644
--- a/src/patch/ghost_zone.hh
+++ b/src/patch/ghost_zone.hh
@@ -1,5 +1,5 @@
 // ghost_zone.hh -- fill in gridfn data in patch ghost zones
-// $Id$
+// $Header$
 //
 // ***** design notes for ghost zones *****
 // ghost_zone - abstract base class to describe ghost zone of patch
diff --git a/src/patch/grid.cc b/src/patch/grid.cc
index 16dc2cd..896ec46 100644
--- a/src/patch/grid.cc
+++ b/src/patch/grid.cc
@@ -1,5 +1,5 @@
 // grid.cc -- classes for a 2D uniform tensor-product grid
-// $Id$
+// $Header$
 //
 // grid_arrays::grid_arrays
 // grid_arrays::setup_gridfn_storage
diff --git a/src/patch/grid.hh b/src/patch/grid.hh
index dfb1364..8da644c 100644
--- a/src/patch/grid.hh
+++ b/src/patch/grid.hh
@@ -1,5 +1,5 @@
 // grid.hh -- classes for a 2D uniform tensor-product grid
-// $Id$
+// $Header$
 //
 // grid_arrays - data arrays for a 2D tensor-product grid
 // grid - uniform 2D tensor-product grid
diff --git a/src/patch/patch.cc b/src/patch/patch.cc
index dde6dcd..b07778b 100644
--- a/src/patch/patch.cc
+++ b/src/patch/patch.cc
@@ -1,5 +1,5 @@
 // patch.cc -- describes a coordinate/grid patch
-// $Id$
+// $Header$
 
 //
 // patch::patch
@@ -8,16 +8,19 @@
 // x_patch::x_patch
 // y_patch::y_patch
 //
+// patch::Jacobian_xyz_wrt_rho_sigma
+//
 // patch::decode_integration_method
 // patch::integrate_gridfn
-// patch::integration_coeff
+/// patch::integration_coeff
+//
 // patch::ghost_zone_on_edge
 // patch::corner_ghost_zone_containing_point
 // patch::ghost_zone_containing_point
 // patch::create_mirror_symmetry_ghost_zone
 // patch::create_periodic_symmetry_ghost_zone
 // patch::create_interpatch_ghost_zone
-// patch::set_ghost_zone
+/// patch::set_ghost_zone
 // patch::edge_adjacent_to_patch
 // patch::assert_all_ghost_zones_fully_setup
 //
@@ -29,6 +32,7 @@
 using std::fprintf;
 using std::printf;
 using std::strcmp;
+using std::sqrt;
 
 #include "cctk.h"
 
@@ -161,6 +165,53 @@ y_patch::y_patch(patch_system &my_patch_system_in, int patch_number_in,
 //******************************************************************************
 
 //
+// This function computes the Jacobian of (x,y,z) with respect to (rho,sigma),
+//	det [g_{ij} (partial x^i/partial y^u) (partial y^j/partial y^v)]
+//
+fp patch::Jacobian_xyz_wrt_rho_sigma(fp r, fp rho, fp sigma,
+				     fp g_xx, fp g_xy, fp g_xz,
+					      fp g_yy, fp g_yz,
+						       fp g_zz)
+	const
+{
+fp partial_x_wrt_rho, partial_x_wrt_sigma;
+fp partial_y_wrt_rho, partial_y_wrt_sigma;
+fp partial_z_wrt_rho, partial_z_wrt_sigma;
+partial_xyz_wrt_rho_sigma(r, rho, sigma,
+			  partial_x_wrt_rho, partial_x_wrt_sigma,
+			  partial_y_wrt_rho, partial_y_wrt_sigma,
+			  partial_z_wrt_rho, partial_z_wrt_sigma);
+
+const fp g_rho_rho = +     partial_x_wrt_rho*partial_x_wrt_rho*g_xx
+		     + 2.0*partial_x_wrt_rho*partial_y_wrt_rho*g_xy
+		     + 2.0*partial_x_wrt_rho*partial_z_wrt_rho*g_xz
+		     +     partial_y_wrt_rho*partial_y_wrt_rho*g_yy
+		     + 2.0*partial_y_wrt_rho*partial_z_wrt_rho*g_yz
+		     +     partial_z_wrt_rho*partial_z_wrt_rho*g_zz;
+const fp g_rho_sigma = +     partial_x_wrt_rho*partial_x_wrt_sigma    *g_xx
+		       + (   partial_x_wrt_rho*partial_y_wrt_sigma
+			   + partial_y_wrt_rho*partial_x_wrt_sigma   )*g_xy
+		       + (   partial_x_wrt_rho*partial_z_wrt_sigma
+			   + partial_z_wrt_rho*partial_x_wrt_sigma   )*g_xz
+		       +     partial_y_wrt_rho*partial_y_wrt_sigma    *g_yy
+		       + (   partial_y_wrt_rho*partial_z_wrt_sigma
+			   + partial_z_wrt_rho*partial_y_wrt_sigma   )*g_yz
+		       +     partial_z_wrt_rho*partial_z_wrt_sigma    *g_zz;
+const fp g_sigma_sigma = +     partial_x_wrt_sigma*partial_x_wrt_sigma*g_xx
+			 + 2.0*partial_x_wrt_sigma*partial_y_wrt_sigma*g_xy
+			 + 2.0*partial_x_wrt_sigma*partial_z_wrt_sigma*g_xz
+			 +     partial_y_wrt_sigma*partial_y_wrt_sigma*g_yy
+			 + 2.0*partial_y_wrt_sigma*partial_z_wrt_sigma*g_yz
+			 +     partial_z_wrt_sigma*partial_z_wrt_sigma*g_zz;
+
+return g_rho_rho*g_sigma_sigma - jtutil::pow2(g_rho_sigma);
+}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
 // This function decodes the character-string name of an integration method
 // into an  enum integration_method .  See the comments in "patch.hh" on the
 // declaration of  enum integration_method  for details on the methods and
@@ -191,18 +242,25 @@ else	error_exit(ERROR_EXIT,
 //
 // This function computes an approximation to the (surface) integral of
 // a gridfn over the patch's nominal area,
-//    area_weighting_flag ? $\int f(\rho,\sigma) \, dA$
-//                        : $\int f(\rho,\sigma) \, d\rho \, d\sigma$
-// where $f(\rho,\sigma)$ is the gridfn and $dA$ is the area element in
-// $(\rho,sigma)$ coordinates.
+//	$\int f(\rho,\sigma) \, dA$
+//		= \int f(\rho,\sigma) \sqrt{|J|} \, d\rho \, d\sigma$
+// where $J$ is the Jacobian of $(x,y,z)$ with respect to $(rho,sigma).
 //
-// The integration scheme is selected based on the  method  argument.
+// Arguments:
+// src_gfn = (in) The gridfn to be integrated.
+// ghosted_radius_gfn = (in) The surface radius.
+// g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
+// method = (in) Selects the integration scheme.
 //
 // Bugs:
-// The way the integration coefficients are computed is somewhat inefficient.
+// - The integration coefficients are computed in a rather inefficient
+//   manner.
 //
 fp patch::integrate_gridfn(int src_gfn,
-			   bool area_weighting_flag,
+			   int ghosted_radius_gfn,
+			   int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					 int g_yy_gfn, int g_yz_gfn,
+						       int g_zz_gfn,
 			   enum integration_method method)
 	const
 {
@@ -211,13 +269,32 @@ fp sum = 0.0;
 	{
 	for (int isigma = min_isigma() ; isigma <= max_isigma() ; ++isigma)
 	{
-	sum += gridfn(src_gfn, irho,isigma)
-	       * integration_coeff(method,
-				   max_irho()-min_irho(),
-				   irho      -min_irho())
-	       * integration_coeff(method,
-				   max_isigma()-min_isigma(),
-				   isigma      -min_isigma());
+	const fp fn = gridfn(src_gfn, irho,isigma);
+
+	const fp r = ghosted_gridfn(ghosted_radius_gfn, irho,isigma);
+	const fp rho   = rho_of_irho    (irho);
+	const fp sigma = sigma_of_isigma(isigma);
+
+	const fp g_xx = gridfn(g_xx_gfn, irho,isigma);
+	const fp g_xy = gridfn(g_xy_gfn, irho,isigma);
+	const fp g_xz = gridfn(g_xz_gfn, irho,isigma);
+	const fp g_yy = gridfn(g_yy_gfn, irho,isigma);
+	const fp g_yz = gridfn(g_yz_gfn, irho,isigma);
+	const fp g_zz = gridfn(g_zz_gfn, irho,isigma);
+
+	const fp Jac = Jacobian_xyz_wrt_rho_sigma(r, rho, sigma,
+						  g_xx, g_xy, g_xz,
+							g_yy, g_yz,
+							      g_zz);
+
+	const fp coeff_rho   = integration_coeff(method,
+						 max_irho()-min_irho(),
+						 irho      -min_irho());
+	const fp coeff_sigma = integration_coeff(method,
+						 max_isigma()-min_isigma(),
+						 isigma      -min_isigma());
+
+	sum += fn * sqrt(jtutil::abs(Jac)) * coeff_rho*coeff_sigma;
 	}
 	}
 return delta_rho() * delta_sigma() * sum;
diff --git a/src/patch/patch.hh b/src/patch/patch.hh
index 025aecd..94b4c0b 100644
--- a/src/patch/patch.hh
+++ b/src/patch/patch.hh
@@ -1,5 +1,5 @@
 // patch.hh -- describes a coordinate/grid patch
-// $Id$
+// $Header$
 //
 // ***** how patch boundaries are handled *****
 // patch - abstract base class to describe a coordinate/grid patch
@@ -218,6 +218,14 @@ public:
 					 fp& xcos, fp& ycos, fp& zcos)
 		const = 0;
 
+	// Jacobian of (x,y,z) with respect to (rho,sigma):
+	// det [g_{ij} (partial x^i/partial y^u) (partial y^j/partial y^v)]
+	fp Jacobian_xyz_wrt_rho_sigma(fp r, fp rho, fp sigma,
+				      fp g_xx, fp g_xy, fp g_xz,
+					       fp g_yy, fp g_yz,
+							fp g_zz)
+		const;
+
 	// partial (x,y,z) / partial (rho,sigma)
 	virtual void partial_xyz_wrt_rho_sigma
 		(fp r, fp rho, fp sigma,
@@ -340,22 +348,23 @@ public:
 	  enum integration_method
 	    decode_integration_method(const char method_string[]);
 
-
-	// integrate a gridfn: computes an approximation to the surface
-	// integral
-	//    area_weighting_flag ? $\int f(\rho,\sigma) \, dA$
-	//                        : $\int f(\rho,\sigma) \, d\rho \, d\sigma$
-	// where $dA$ is the area element in $(\rho,sigma)$ coordinates
+	// compute the surface integral of a gridfn over the patch,
+	//	$\int f(\rho,\sigma) \, dA$
+	//		= \int f(\rho,\sigma) \sqrt{|J|} \, d\rho \, d\sigma$
+	// where $J$ is the Jacobian of $(x,y,z)$ with respect to $(rho,sigma)
 	// ... integration method selected by  method  argument
-	//     FIXME: right now this is ignored :( :(
 	fp integrate_gridfn(int src_gfn,
-			    bool area_weighting_flag,
+			    int ghosted_radius_gfn,
+			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					  int g_yy_gfn, int g_yz_gfn,
+							int g_zz_gfn,
 			    enum integration_method method)
 		const;
-private:
+
 	// compute integration coefficient $c_i$ where
 	// $\int_{x_0}^{x_N} f(x) \, dx
 	//	\approx \Delta x \, \sum_{i=0}^N c_i f(x_i)$
+private:
 	static
 	  fp integration_coeff(enum integration_method method, int N, int i);
 
diff --git a/src/patch/patch_edge.hh b/src/patch/patch_edge.hh
index b1586b5..4fba1b9 100644
--- a/src/patch/patch_edge.hh
+++ b/src/patch/patch_edge.hh
@@ -1,5 +1,5 @@
 // patch_edge.hh -- perpendicular/parallel geometry of one side of a patch
-// $Id$
+// $Header$
 //
 
 //
diff --git a/src/patch/patch_info.cc b/src/patch/patch_info.cc
index 7d6488e..db49ba7 100644
--- a/src/patch/patch_info.cc
+++ b/src/patch/patch_info.cc
@@ -1,5 +1,5 @@
 // patch_info.cc -- POD struct of minimal info varying from one patch to another
-// $Id$
+// $Header$
 
 //
 // patch_info::grid_array_pars
diff --git a/src/patch/patch_info.hh b/src/patch/patch_info.hh
index 7f58f21..f1e0043 100644
--- a/src/patch/patch_info.hh
+++ b/src/patch/patch_info.hh
@@ -1,5 +1,5 @@
 // patch_info.hh -- POD struct of minimal info varying from one patch to another
-// $Id$
+// $Header$
 
 //
 // prerequisites:
diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index 2005d89..a0ad72f 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -1,5 +1,5 @@
 // patch_interp.cc -- interpolation from a patch
-// $Id$
+// $Header$
 
 //
 // patch_interp::patch_interp
diff --git a/src/patch/patch_interp.hh b/src/patch/patch_interp.hh
index 8ad3357..b4bff04 100644
--- a/src/patch/patch_interp.hh
+++ b/src/patch/patch_interp.hh
@@ -1,5 +1,5 @@
 // patch_interp.hh -- interpolation from a patch
-// $Id$
+// $Header$
 
 //
 // prerequisites:
diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index c8ed7d0..ee80535 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -1,5 +1,5 @@
 // patch_system.cc -- describes the (an) entire system of patches
-// $Id$
+// $Header$
 //
 // patch_system::patch_system
 // patch_system::~patch_system
@@ -1257,19 +1257,22 @@ norms.reset();
 
 //
 // This function computes an approximation to the (surface) integral of
-// a gridfn over the whole patch system,
-//    area_weighting_flag ? $\int f(\rho,\sigma) \, dA$
-//                        : $\int f(\rho,\sigma) \, d\rho \, d\sigma$
-// where $f(\rho,\sigma)$ is the gridfn and $dA$ is the area element in
-// $(\rho,sigma)$ coordinates.
+// a gridfn over the patch system,
+//	$\int f(\rho,\sigma) \, dA$
+//		= \int f(\rho,\sigma) \sqrt{|J|} \, d\rho \, d\sigma$
+// where $J$ is the Jacobian of $(x,y,z)$ with respect to $(rho,sigma).
 //
-// The integration scheme is selected based on the  method  argument.
-//
-// Bugs:
-// The way the integration coefficients are computed is somewhat inefficient.
+// Arguments:
+// src_gfn = (in) The gridfn to be integrated.
+// ghosted_radius_gfn = (in) The surface radius.
+// g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
+// method = (in) Selects the integration scheme.
 //
 fp patch_system::integrate_gridfn(int src_gfn,
-				  bool area_weighting_flag,
+				  int ghosted_radius_gfn,
+				  int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+						int g_yy_gfn, int g_yz_gfn,
+							      int g_zz_gfn,
 				  enum patch::integration_method method)
 	const
 {
@@ -1277,7 +1280,21 @@ fp integral = 0.0;
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
 	const patch& p = ith_patch(pn);
-	integral += p.integrate_gridfn(src_gfn, area_weighting_flag, method);
+	integral += p.integrate_gridfn(src_gfn,
+				       ghosted_radius_gfn,
+				       g_xx_gfn, g_xy_gfn, g_xz_gfn,
+						 g_yy_gfn, g_yz_gfn,
+							   g_zz_gfn,
+				       method);
+/* PICKLE */
+	const fp area = p.integrate_gridfn(src_gfn,
+				       ghosted_radius_gfn,
+				       g_xx_gfn, g_xy_gfn, g_xz_gfn,
+						 g_yy_gfn, g_yz_gfn,
+							   g_zz_gfn,
+				       method);
+	if (src_gfn == max_gfn())
+	   then printf("%s patch: area=%.10f\n", p.name(), area);
 	}
 return integral;
 }
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index 75ca8cc..2b2f1c9 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -1,5 +1,5 @@
 // patch_system.hh -- describes the (an) entire system of interlinked patches
-// $Id$
+// $Header$
 //
 // patch_system - describes a system of interlinked patches
 //
@@ -260,14 +260,16 @@ public:
 	void ghosted_gridfn_norms(int ghosted_src_gfn, jtutil::norm<fp>& norms)
 		const;
 
-	// integrate a gridfn: computes an approximation to the surface
-	// integral over the whole patch system
-	//    area_weighting_flag ? $\int f(\rho,\sigma) \, dA$
-	//                        : $\int f(\rho,\sigma) \, d\rho \, d\sigma$
-	// where $dA$ is the area element in $(\rho,sigma)$ coordinates
+	// compute the surface integral of a gridfn over the patch system,
+	//	$\int f(\rho,\sigma) \, dA$
+	//		= \int f(\rho,\sigma) \sqrt{|J|} \, d\rho \, d\sigma$
+	// where $J$ is the Jacobian of $(x,y,z)$ with respect to $(rho,sigma)
 	// ... integration method selected by  method  argument
 	fp integrate_gridfn(int src_gfn,
-			    bool area_weighting_flag,
+			    int ghosted_radius_gfn,
+			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					  int g_yy_gfn, int g_yz_gfn,
+							int g_zz_gfn,
 			    enum patch::integration_method method)
 		const;
 
diff --git a/src/patch/patch_system_info.hh b/src/patch/patch_system_info.hh
index dd6cf9f..79f2b7a 100644
--- a/src/patch/patch_system_info.hh
+++ b/src/patch/patch_system_info.hh
@@ -1,5 +1,5 @@
 // patch_system_info.hh -- static info describing various types of patch systems
-// $Id$
+// $Header$
 
 //
 // prerequisites:
diff --git a/src/patch/test_coords.cc b/src/patch/test_coords.cc
index 2fc47de..05e3628 100644
--- a/src/patch/test_coords.cc
+++ b/src/patch/test_coords.cc
@@ -1,5 +1,5 @@
 // test_coords.cc -- test driver for coordinate systems/conversions
-// $Id$
+// $Header$
 
 #include <stdio.h>
 #include <assert.h>
diff --git a/src/patch/test_coords2.cc b/src/patch/test_coords2.cc
index 0dcf056..762ad14 100644
--- a/src/patch/test_coords2.cc
+++ b/src/patch/test_coords2.cc
@@ -1,5 +1,5 @@
 // test_coords2.cc -- test driver #2 for coordinate systems/conversions
-// $Id$
+// $Header$
 
 #include <stdio.h>
 #include <assert.h>
diff --git a/src/patch/test_patch_system.cc b/src/patch/test_patch_system.cc
index d4579df..d4610f2 100644
--- a/src/patch/test_patch_system.cc
+++ b/src/patch/test_patch_system.cc
@@ -1,5 +1,5 @@
 // test_patch_system.cc -- test driver for patch_system::
-// $Id$
+// $Header$
 
 //
 // <<<misc global data>>>
diff --git a/src/patch/test_patch_system.maple b/src/patch/test_patch_system.maple
index 17cacef..3315922 100644
--- a/src/patch/test_patch_system.maple
+++ b/src/patch/test_patch_system.maple
@@ -1,5 +1,5 @@
 # test_patch_system.maple -- test function for  test_fd_grid
-# $Id$
+# $Header$
 
 # test function
 fn := exp(sin(1.38*rho)) * tanh(0.17+0.83*sin(sigma)^2);