add code to compute surface integrals of gridfns over patches

and over the whole patch system
-- note the volume element isn't included yet


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

1 files changed, 101 insertions, 0 deletions

diff --git a/src/patch/patch.hh b/src/patch/patch.hh
index 8a90d56..950a7f6 100644
--- a/src/patch/patch.hh
+++ b/src/patch/patch.hh
@@ -252,6 +252,107 @@ public:
 
 
 	//
+	// ***** gridfn operations *****
+	//
+public:
+
+	//
+	// The following enum describes the integration methods supported
+	// by  integrate_gridfn() .
+	//
+	// For convenience of exposition we describe the methods as if for
+	// 1-D integration, but  integrate_gridfn()  actually does 2-D
+	// (surface) integration over the patch.
+	//
+	// Suppose we're computing $\int_{x_0}^{x^N} f(x) \, dx$, using the
+	// equally spaced integration points $f_0$, $f_1$, \dots, $f_N$,
+	// spaced $\Delta x$ apart.  Then the integration methods are as
+	// follows, with the convention that $\langle X \rangle$ denotes
+	// indefinite repetition of the "X" terms, depending on N:
+	//
+	enum	integration_method
+		{
+		// Trapezoid rule
+		// ... character-string name "trapezoid" or "trapezoid rule"
+		// ... 2nd order accurate for smooth functions
+		// ... requires N >= 1
+		// $$
+		// \Delta x \left[
+		//            \half f_0
+		//          + \langle
+		//            f_k
+		//            \rangle
+		//          + \half f_N
+		//          \right]
+		// $$
+		integration_method__trapezoid,
+
+		// Simpson's rule
+		// ... character-string name "Simpson" or "Simpson's rule"
+		// ... 4th order accurate for smooth functions
+		// ... requires N >= 2 and N even
+		// $$
+		// \Delta x \left[
+		//            \frac{1}{3} f_0
+		//          + \frac{4}{3} f_1
+		//          + \langle
+		//            \frac{2}{3} f_{2k} + \frac{4}{3} f_{2k+1}
+		//            \rangle
+		//          + \frac{1}{3} f_N
+		//          \right]
+		// $$
+		integration_method__Simpson,
+
+		// Simpson's rule, variant form
+		// ... characgter-string name "Simpson (variant)"
+		//     or "Simpson's rule (variant)"
+		// ... described in Numerical Recipes 1st edition (4.1.14)
+		// ... 4th order accurate for smooth functions
+		// ... requires N >= 7
+		// $$
+		// \Delta x \left[
+		//            \frac{17}{48} f_0
+		//          + \frac{59}{48} f_1
+		//          + \frac{43}{48} f_2
+		//          + \frac{49}{48} f_3
+		//          + \langle
+		//            f_k
+		//            \rangle
+		//          + \frac{49}{48} f_{N-3}
+		//          + \frac{43}{48} f_{N-2}
+		//          + \frac{59}{48} f_{N-1}
+		//          + \frac{17}{48} f_N
+		//          \right]
+		// $$
+		integration_method__Simpson_variant // no comma here!
+		};
+
+	// decode character string name into internal enum
+	static
+	  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
+	// ... integration method selected by  method  argument
+	//     FIXME: right now this is ignored :( :(
+	fp integrate_gridfn(int src_gfn,
+			    bool area_weighting_flag,
+			    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)$
+	static
+	  fp integration_coeff(enum integration_method method, int N, int i);
+
+
+	//
 	// ***** patch edges ****
 	//
 public: