add support for computing line integrals around a surface

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

2 files changed, 235 insertions, 48 deletions

diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index 0aefef6..1b3ee83 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -34,6 +34,11 @@
 // patch_system::add_to_ghosted_gridfn
 // patch_system::gridfn_norms
 // patch_system::ghosted_gridfn_norms
+//
+// patch_system::xy_circumference
+// patch_system::xz_circumference
+// patch_system::yz_circumference
+// patch_system::arc_length_sum
 // patch_system::integrate_gridfn
 //
 // patch_system::print_unknown_gridfn
@@ -163,8 +168,7 @@ if (print_summary_msg_flag)
 switch	(type_in)
 	{
 case patch_system__full_sphere:
-	create_patches(patch_system_info::full_sphere
-					::patch_info_array,
+	create_patches(patch_system_info::full_sphere::patch_info_array,
 		       ghost_zone_width, patch_extend_width,
 		       N_zones_per_right_angle,
 		       print_detailed_msg_flag);
@@ -178,8 +182,7 @@ case patch_system__full_sphere:
 	break;
 
 case patch_system__plus_z_hemisphere:
-	create_patches(patch_system_info::plus_z_hemisphere
-					::patch_info_array,
+	create_patches(patch_system_info::plus_z_hemisphere::patch_info_array,
 		       ghost_zone_width, patch_extend_width,
 		       N_zones_per_right_angle,
 		       print_detailed_msg_flag);
@@ -193,8 +196,7 @@ case patch_system__plus_z_hemisphere:
 	break;
 
 case patch_system__plus_xy_quadrant_mirrored:
-	create_patches(patch_system_info::plus_xy_quadrant_mirrored
-					::patch_info_array,
+	create_patches(patch_system_info::plus_xy_quadrant::patch_info_array,
 		       ghost_zone_width, patch_extend_width,
 		       N_zones_per_right_angle,
 		       print_detailed_msg_flag);
@@ -208,8 +210,7 @@ case patch_system__plus_xy_quadrant_mirrored:
 	break;
 
 case patch_system__plus_xy_quadrant_rotating:
-	create_patches(patch_system_info::plus_xy_quadrant_rotating
-					::patch_info_array,
+	create_patches(patch_system_info::plus_xy_quadrant::patch_info_array,
 		       ghost_zone_width, patch_extend_width,
 		       N_zones_per_right_angle,
 		       print_detailed_msg_flag);
@@ -223,8 +224,7 @@ case patch_system__plus_xy_quadrant_rotating:
 	break;
 
 case patch_system__plus_xz_quadrant_rotating:
-	create_patches(patch_system_info::plus_xz_quadrant_rotating
-					::patch_info_array,
+	create_patches(patch_system_info::plus_xz_quadrant::patch_info_array,
 		       ghost_zone_width, patch_extend_width,
 		       N_zones_per_right_angle,
 		       print_detailed_msg_flag);
@@ -238,8 +238,7 @@ case patch_system__plus_xz_quadrant_rotating:
 	break;
 
 case patch_system__plus_xyz_octant_mirrored:
-	create_patches(patch_system_info::plus_xyz_octant_mirrored
-					::patch_info_array,
+	create_patches(patch_system_info::plus_xyz_octant::patch_info_array,
 		       ghost_zone_width, patch_extend_width,
 		       N_zones_per_right_angle,
 		       print_detailed_msg_flag);
@@ -253,8 +252,7 @@ case patch_system__plus_xyz_octant_mirrored:
 	break;
 
 case patch_system__plus_xyz_octant_rotating:
-	create_patches(patch_system_info::plus_xyz_octant_rotating
-					::patch_info_array,
+	create_patches(patch_system_info::plus_xyz_octant::patch_info_array,
 		       ghost_zone_width, patch_extend_width,
 		       N_zones_per_right_angle,
 		       print_detailed_msg_flag);
@@ -624,12 +622,12 @@ if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "      seting up full sphere ghost zones");
 
-patch& pz = ith_patch(patch_number_of_name("+z"));
-patch& px = ith_patch(patch_number_of_name("+x"));
-patch& py = ith_patch(patch_number_of_name("+y"));
-patch& mx = ith_patch(patch_number_of_name("-x"));
-patch& my = ith_patch(patch_number_of_name("-y"));
-patch& mz = ith_patch(patch_number_of_name("-z"));
+patch& pz = ith_patch(patch_system_info::full_sphere::patch_number__pz);
+patch& px = ith_patch(patch_system_info::full_sphere::patch_number__px);
+patch& py = ith_patch(patch_system_info::full_sphere::patch_number__py);
+patch& mx = ith_patch(patch_system_info::full_sphere::patch_number__mx);
+patch& my = ith_patch(patch_system_info::full_sphere::patch_number__my);
+patch& mz = ith_patch(patch_system_info::full_sphere::patch_number__mz);
 
 // create the ghost zones
 if (print_msg_flag)
@@ -707,11 +705,11 @@ if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "      setting up +z hemisphere ghost zones");
 
-patch& pz = ith_patch(patch_number_of_name("+z"));
-patch& px = ith_patch(patch_number_of_name("+x"));
-patch& py = ith_patch(patch_number_of_name("+y"));
-patch& mx = ith_patch(patch_number_of_name("-x"));
-patch& my = ith_patch(patch_number_of_name("-y"));
+patch& pz = ith_patch(patch_system_info::plus_z_hemisphere::patch_number__pz);
+patch& px = ith_patch(patch_system_info::plus_z_hemisphere::patch_number__px);
+patch& py = ith_patch(patch_system_info::plus_z_hemisphere::patch_number__py);
+patch& mx = ith_patch(patch_system_info::plus_z_hemisphere::patch_number__mx);
+patch& my = ith_patch(patch_system_info::plus_z_hemisphere::patch_number__my);
 
 // create the ghost zones
 if (print_msg_flag)
@@ -777,10 +775,10 @@ if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "      setting up +xy quadrant (mirrored) ghost zones");
 
-patch& pz = ith_patch(patch_number_of_name("+z"));
-patch& px = ith_patch(patch_number_of_name("+x"));
-patch& py = ith_patch(patch_number_of_name("+y"));
-patch& mz = ith_patch(patch_number_of_name("-z"));
+patch& pz = ith_patch(patch_system_info::plus_xy_quadrant::patch_number__pz);
+patch& px = ith_patch(patch_system_info::plus_xy_quadrant::patch_number__px);
+patch& py = ith_patch(patch_system_info::plus_xy_quadrant::patch_number__py);
+patch& mz = ith_patch(patch_system_info::plus_xy_quadrant::patch_number__mz);
 
 // create the ghost zones
 if (print_msg_flag)
@@ -836,10 +834,10 @@ if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "      setting up +xy quadrant (rotating) ghost zones");
 
-patch& pz = ith_patch(patch_number_of_name("+z"));
-patch& px = ith_patch(patch_number_of_name("+x"));
-patch& py = ith_patch(patch_number_of_name("+y"));
-patch& mz = ith_patch(patch_number_of_name("-z"));
+patch& pz = ith_patch(patch_system_info::plus_xy_quadrant::patch_number__pz);
+patch& px = ith_patch(patch_system_info::plus_xy_quadrant::patch_number__px);
+patch& py = ith_patch(patch_system_info::plus_xy_quadrant::patch_number__py);
+patch& mz = ith_patch(patch_system_info::plus_xy_quadrant::patch_number__mz);
 
 // create the ghost zones
 if (print_msg_flag)
@@ -898,10 +896,10 @@ if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "      setting up +xz quadrant (rotating) ghost zones");
 
-patch& pz = ith_patch(patch_number_of_name("+z"));
-patch& px = ith_patch(patch_number_of_name("+x"));
-patch& py = ith_patch(patch_number_of_name("+y"));
-patch& my = ith_patch(patch_number_of_name("-y"));
+patch& pz = ith_patch(patch_system_info::plus_xz_quadrant::patch_number__pz);
+patch& px = ith_patch(patch_system_info::plus_xz_quadrant::patch_number__px);
+patch& py = ith_patch(patch_system_info::plus_xz_quadrant::patch_number__py);
+patch& my = ith_patch(patch_system_info::plus_xz_quadrant::patch_number__my);
 
 // create the ghost zones
 if (print_msg_flag)
@@ -960,9 +958,9 @@ if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "      setting up +xyz octant (mirrored) ghost zones");
 
-patch& pz = ith_patch(patch_number_of_name("+z"));
-patch& px = ith_patch(patch_number_of_name("+x"));
-patch& py = ith_patch(patch_number_of_name("+y"));
+patch& pz = ith_patch(patch_system_info::plus_xyz_octant::patch_number__pz);
+patch& px = ith_patch(patch_system_info::plus_xyz_octant::patch_number__px);
+patch& py = ith_patch(patch_system_info::plus_xyz_octant::patch_number__py);
 
 // create the ghost zones
 if (print_msg_flag)
@@ -1010,9 +1008,9 @@ if (print_msg_flag)
    then CCTK_VInfo(CCTK_THORNSTRING,
 		   "      setting up +xyz octant (rotating) ghost zones");
 
-patch& pz = ith_patch(patch_number_of_name("+z"));
-patch& px = ith_patch(patch_number_of_name("+x"));
-patch& py = ith_patch(patch_number_of_name("+y"));
+patch& pz = ith_patch(patch_system_info::plus_xyz_octant::patch_number__pz);
+patch& px = ith_patch(patch_system_info::plus_xyz_octant::patch_number__px);
+patch& py = ith_patch(patch_system_info::plus_xyz_octant::patch_number__py);
 
 // create the ghost zones
 if (print_msg_flag)
@@ -1143,15 +1141,13 @@ case patch_system__full_sphere:
 case patch_system__plus_z_hemisphere:
 	return patch_system_info::plus_z_hemisphere::N_patches;
 case patch_system__plus_xy_quadrant_mirrored:
-	return patch_system_info::plus_xy_quadrant_mirrored::N_patches;
 case patch_system__plus_xy_quadrant_rotating:
-	return patch_system_info::plus_xy_quadrant_rotating::N_patches;
+	return patch_system_info::plus_xy_quadrant::N_patches;
 case patch_system__plus_xz_quadrant_rotating:
-	return patch_system_info::plus_xz_quadrant_rotating::N_patches;
+	return patch_system_info::plus_xz_quadrant::N_patches;
 case patch_system__plus_xyz_octant_mirrored:
-	return patch_system_info::plus_xyz_octant_mirrored::N_patches;
 case patch_system__plus_xyz_octant_rotating:
-	return patch_system_info::plus_xyz_octant_rotating::N_patches;
+	return patch_system_info::plus_xyz_octant::N_patches;
 default:
 	error_exit(PANIC_EXIT,
 "***** patch_system::N_patches_of_type(): bad type=(int)%d!\n"
@@ -1543,6 +1539,170 @@ norms.reset();
 }
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function computes an approximation to the circumference of a
+// surface in the xy plane.  That is, we integrate in the phi direction
+// around the equator (the +/- x and +/- y patches).  Note that we
+// copute the full circumference all around the sphere, even if the
+// patch system only covers a proper subset of this.
+//
+// Arguments:
+// 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::xy_circumference(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
+{
+fp circumference = 0.0;
+
+    for (int pn = 0 ; pn < N_patches() ; ++pn)
+    {
+    const patch& p = ith_patch(pn);
+    if ((p.ctype() == 'x') || (p.ctype() == 'y'))
+       then circumference += p.plane_arc_length('z',
+						ghosted_radius_gfn,
+						g_xx_gfn, g_xy_gfn, g_xz_gfn,
+							  g_yy_gfn, g_yz_gfn,
+								    g_zz_gfn,
+						method);
+    }
+
+switch	(type())
+	{
+case patch_system__full_sphere:						break;
+case patch_system__plus_z_hemisphere:					break;
+case patch_system__plus_xy_quadrant_mirrored:	circumference *= 4.0;	break;
+case patch_system__plus_xy_quadrant_rotating:	circumference *= 4.0;	break;
+case patch_system__plus_xz_quadrant_rotating:	circumference *= 2.0;	break;
+case patch_system__plus_xyz_octant_mirrored:	circumference *= 4.0;	break;
+case patch_system__plus_xyz_octant_rotating:	circumference *= 4.0;	break;
+default:
+	error_exit(PANIC_EXIT,
+"***** patch_system::xy_circumference(): unknown type=(int)%d!\n"
+"                                        (this should never happen!)\n",
+		   int(type()));				/*NOTREACHED*/
+	}
+
+return circumference;
+}
+
+//******************************************************************************
+
+//
+// This function coputes an approximation to the circumference of a
+// surface in the xz plane.  That is, we integrate in the nu direction
+// around the y=0 plane (the +/- x and +/- z patches).  Note that we
+// copute the full circumference all around the sphere, even if the
+// patch system only covers a proper subset of this.
+//
+// Arguments:
+// 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::xz_circumference(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
+{
+fp circumference = 0.0;
+
+    for (int pn = 0 ; pn < N_patches() ; ++pn)
+    {
+    const patch& p = ith_patch(pn);
+    if ((p.ctype() == 'x') || (p.ctype() == 'z'))
+       then circumference += p.plane_arc_length('y',
+						ghosted_radius_gfn,
+						g_xx_gfn, g_xy_gfn, g_xz_gfn,
+							  g_yy_gfn, g_yz_gfn,
+								    g_zz_gfn,
+						method);
+    }
+
+switch	(type())
+	{
+case patch_system__full_sphere:						break;
+case patch_system__plus_z_hemisphere:		circumference *= 2.0;	break;
+case patch_system__plus_xy_quadrant_mirrored:	circumference *= 2.0;	break;
+case patch_system__plus_xy_quadrant_rotating:	circumference *= 2.0;	break;
+case patch_system__plus_xz_quadrant_rotating:	circumference *= 4.0;	break;
+case patch_system__plus_xyz_octant_mirrored:	circumference *= 4.0;	break;
+case patch_system__plus_xyz_octant_rotating:	circumference *= 4.0;	break;
+default:
+	error_exit(PANIC_EXIT,
+"***** patch_system::xz_circumference(): unknown type=(int)%d!\n"
+"                                        (this should never happen!)\n",
+		   int(type()));				/*NOTREACHED*/
+	}
+
+return circumference;
+}
+
+//******************************************************************************
+
+//
+// This function coputes an approximation to the circumference of a
+// surface in the yz plane.  That is, we integrate in the mu direction
+// around the x=0 plane (the +/- y and +/- z patches).  Note that we
+// copute the full circumference all around the sphere, even if the
+// patch system only covers a proper subset of this.
+//
+// Arguments:
+// ghosted_radius_gfn = (in) The surface radius.
+// g_{xx,xy,yz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
+// method = (in) Selects the integration scheme.
+//
+fp patch_system::yz_circumference(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
+{
+fp circumference = 0.0;
+
+    for (int pn = 0 ; pn < N_patches() ; ++pn)
+    {
+    const patch& p = ith_patch(pn);
+    if ((p.ctype() == 'y') || (p.ctype() == 'z'))
+       then circumference += p.plane_arc_length('x',
+						ghosted_radius_gfn,
+						g_xx_gfn, g_xy_gfn, g_xz_gfn,
+							  g_yy_gfn, g_yz_gfn,
+								    g_zz_gfn,
+						method);
+    }
+
+switch	(type())
+	{
+case patch_system__full_sphere:						break;
+case patch_system__plus_z_hemisphere:		circumference *= 2.0;	break;
+case patch_system__plus_xy_quadrant_mirrored:	circumference *= 2.0;	break;
+case patch_system__plus_xy_quadrant_rotating:	circumference *= 2.0;	break;
+case patch_system__plus_xz_quadrant_rotating:	circumference *= 2.0;	break;
+case patch_system__plus_xyz_octant_mirrored:	circumference *= 4.0;	break;
+case patch_system__plus_xyz_octant_rotating:	circumference *= 4.0;	break;
+default:
+	error_exit(PANIC_EXIT,
+"***** patch_system::xz_circumference(): unknown type=(int)%d!\n"
+"                                        (this should never happen!)\n",
+		   int(type()));				/*NOTREACHED*/
+	}
+
+return circumference;
+}
+
+//******************************************************************************
 
 //
 // This function computes an approximation to the (surface) integral of
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index d4fe3ce..24e61de 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -278,6 +278,33 @@ public:
 	void ghosted_gridfn_norms(int ghosted_src_gfn, jtutil::norm<fp>& norms)
 		const;
 
+
+	//
+	// ***** line/surface integrals *****
+	//
+
+	// compute the circumference of a surface in the {xy, xz, yz} plane
+	// ... note this is the full circumference all around the sphere,
+	//     even if the patch system only covers a proper subset of this
+	fp xy_circumference(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;
+	fp xz_circumference(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;
+	fp yz_circumference(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;
+
 	// 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$