add explicit gridfns for partial_k g_ij

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

4 files changed, 28 insertions, 11 deletions

diff --git a/src/gr/gr_gfas.minc b/src/gr/gr_gfas.minc
index cb253f4..df3f578 100644
--- a/src/gr/gr_gfas.minc
+++ b/src/gr/gr_gfas.minc
@@ -8,6 +8,8 @@ K_uu, K_uu__fnd,
 
 K, K__fnd,
 
+partial_d_g_dd,
+
 partial_d_ln_sqrt_g, partial_d_ln_sqrt_g__fnd,
 partial_d_g_uu, partial_d_g_uu__fnd,
 
diff --git a/src/gr/metric_derivs.maple b/src/gr/metric_derivs.maple
index 4c922f3..2f30e91 100644
--- a/src/gr/metric_derivs.maple
+++ b/src/gr/metric_derivs.maple
@@ -16,11 +16,14 @@
 #
 # Inputs:
 #	N
-#	g_dd
+#	partial_d_g_dd
 #	g_uu
 #
 # Outputs (Maple + C code file "metric_det_gradient.c")
-#	partial_d_ln_sqrt_g = partial_d_ln_sqrt_g__fnd(g_uu, g_dd)
+#	partial_d_ln_sqrt_g = partial_d_ln_sqrt_g__fnd(g_uu, partial_d_g_dd)
+#
+# Outputs (Maple + C code file "inverse_metric_gradient.c")
+#	partial_d_g_uu = partial_d_g_uu__fnd(g_uu, partial_d_g_dd)
 #
 metric_derivs :=
 proc(cg_flag::boolean)
@@ -40,11 +43,11 @@ end proc;
 #
 # Inputs:
 #	N
-#	g_dd
+#	partial_d_g_dd
 #	g_uu
 #
 # Outputs (Maple + C code file "metric_det_gradient.c")
-#	partial_d_ln_sqrt_g = partial_d_ln_sqrt_g__fnd(g_uu, g_dd)
+#	partial_d_ln_sqrt_g = partial_d_ln_sqrt_g__fnd(g_uu, partial_d_g_dd)
 #
 metric_det_gradient :=
 proc(cg_flag::boolean)
@@ -64,7 +67,7 @@ assert_fnd_exists(partial_d_ln_sqrt_g, fnd);
 	do
 	partial_d_ln_sqrt_g__fnd[k]
 		:= simplify(
-		      (1/2) * msum('g_uu[i,j] * Diff(g_dd[i,j], x_xyz[k])',
+		      (1/2) * msum('g_uu[i,j] * partial_d_g_dd[k,i,j]',
 				   'i'=1..N, 'j'=1..N)
 			   );
 	end do;
@@ -86,11 +89,11 @@ end proc;
 #
 # Inputs:
 #	N
-#	g_dd
+#	partial_d_g_dd
 #	g_uu
 #
 # Outputs (Maple + C code file "inverse_metric_gradient.c")
-#	partial_d_g_uu = partial_d_g_uu__fnd(g_uu, g_dd)
+#	partial_d_g_uu = partial_d_g_uu__fnd(g_uu, partial_d_g_dd)
 #
 inverse_metric_gradient :=
 proc(cg_flag::boolean)
@@ -114,7 +117,7 @@ assert_fnd_exists(partial_d_g_uu, fnd);
 	do
 	partial_d_g_uu__fnd[k,i,j]
 	    := simplify(
-		  - msum('g_uu[i,m] * g_uu[j,n] * Diff(g_dd[m,n], x_xyz[k])',
+		  - msum('g_uu[i,m] * g_uu[j,n] * partial_d_g_dd[k,m,n]',
 			 'm'=1..N, 'n'=1..N)
 		       );
 	end do;
diff --git a/src/gr/setup_gfas.maple b/src/gr/setup_gfas.maple
index 93c74ae..5283d48 100644
--- a/src/gr/setup_gfas.maple
+++ b/src/gr/setup_gfas.maple
@@ -23,11 +23,16 @@ make_gfa('g_uu', {inert,fnd}, [1..N, 1..N], symmetric);
 make_gfa('K_uu', {inert,fnd}, [1..N, 1..N], symmetric);
 make_gfa('K', {inert, fnd}, [], none);
 
-# partial derivatives of metric determinant
+# xyz partial derivatives of metric
+# ... as far as Maple is concerned, these are indeed separate gridfns;
+#     the derivatives are actually taken by the interpolator
+make_gfa('partial_d_g_dd', {inert}, [1..N, 1..N, 1..N], symmetric3_23);
+
+# xyz partial derivatives of metric determinant
 make_gfa('partial_d_ln_sqrt_g', {inert,fnd}, [1..N], none);
 
 # xyz partial derivatives of inverse metric
-# ... actually computed via xyz derivatives of g_dd
+# (computed in terms of xyz derivatives of g_dd)
 make_gfa('partial_d_g_uu', {inert,fnd}, [1..N, 1..N, 1..N], symmetric3_23);
 
 # radius of horizon
diff --git a/src/maple/Diff.maple b/src/maple/Diff.maple
index 7eaf030..3ffec13 100644
--- a/src/maple/Diff.maple
+++ b/src/maple/Diff.maple
@@ -234,6 +234,7 @@ end proc;
 #
 # It currently knows about the following simplifications:
 # - Diff(X_ud[u,i], x_xyz[j]) --> X_udd[u,i,j]
+# - Diff(g_dd[i,j], x_xyz[k]) --> partial_d_g_dd[k,i,j]
 #
 # Anything else is returned unchanged.  (To avoid infinite recursion,
 # such a return is *unevaluated*.)
@@ -248,7 +249,8 @@ proc(operand)			# varargs
 option remember;		# performance optimization
 global
   @include "../maple/coords.minc",
-  @include "../maple/gfa.minc";
+  @include "../maple/gfa.minc",
+  @include "../gr/gr_gfas.minc";
 local var_list, posn;
 
 var_list := [args[2..nargs]];
@@ -258,6 +260,11 @@ if ( type(operand, indexed) and (op(0,operand) = 'X_ud')
    then return X_udd[op(operand), posn];
 end if;
 
+if ( type(operand, indexed) and (op(0,operand) = 'g_dd')
+     and (nops(var_list) = 1) and member(var_list[1],x_xyz_list,'posn') )
+   then return partial_d_g_dd[op(operand), posn];
+end if;
+
 # unevaluated return to avoid infinite recursion
 return 'Diff'(operand, op(var_list));
 end proc;