# Maple code to set up all gridfn arrays
# $Id$

#
# setup_gr_gfas - setup all the GR gfas
# `Diff/gridfn2` - Diff() simplification based on gridfn properties known here
#

################################################################################

setup_gr_gfas :=
proc()
global
  @include "../maple/coords.minc",
  @include "../maple/gfa.minc",
  @include "../gr/gr_gfas.minc";

# metric and extrinsic curvature
make_gfa('g_dd', {inert}, [1..N, 1..N], symmetric);
make_gfa('K_dd', {inert}, [1..N, 1..N], symmetric);

# index-raised and contracted metric and extrinsic curvature
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);

# 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
# (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
make_gfa('h', {inert,fnd}, [], none);

# outward-pointing *non*-unit normal (covector) to horizon
# and it's xyz-coordinate partial derivatives
make_gfa('s_d', {inert,fnd}, [1..N], none);
make_gfa('partial_d_s_d', {inert,fnd}, [1..N, 1..N], none);

# subexpressions for computing LHS of apparent horizon equation
# ... these are defined by (15) in my 1996 apparent horizon finding paper
make_gfa('HA', {inert,fnd}, [], none);
make_gfa('HB', {inert,fnd}, [], none);
make_gfa('HC', {inert,fnd}, [], none);
make_gfa('HD', {inert,fnd}, [], none);

# LHS of apparent horizon equation
make_gfa('H', {inert,fnd}, [], none);

# Jacobian coefficients for H[ABCD]
# these are the partial derivatives of H[ABCD]
# ... wrt Diff(h,x_rs[x])
make_gfa('partial_HA_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
make_gfa('partial_HB_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
make_gfa('partial_HC_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
make_gfa('partial_HD_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
# ... wrt Diff(h,x_rs[x],x_rs[y])
make_gfa('partial_HA_wrt_partial_dd_h', {inert,fnd},
	 [1..N_ang, 1..N_ang], symmetric);
make_gfa('partial_HB_wrt_partial_dd_h', {inert,fnd},
	 [1..N_ang, 1..N_ang], symmetric);

# Jacobian coefficients for H itself
make_gfa('partial_H_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
make_gfa('partial_H_wrt_partial_dd_h', {inert,fnd},
	 [1..N_ang, 1..N_ang], symmetric);

NULL;
end proc;


################################################################################

#
# This function implements further simplification rules for  Diff()
# based on gridfn properties which are known here but not in ../maple/.
#
# It currently knows about the following simplifications:
# - 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*.)
#
# Arguments:
# operand = (in) The thing to be differentiated.
# var_seq = (in) (varargs) An expression sequence of the variables to
#			   differentiate with respect to.
#
`Diff/gridfn2` :=
proc(operand)			# varargs
option remember;		# performance optimization
global
  @include "../maple/coords.minc",
  @include "../maple/gfa.minc",
  @include "../gr/gr_gfas.minc";
local var_list, posn;

var_list := [args[2..nargs]];

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[posn, op(operand)];
end if;

# unevaluated return to avoid infinite recursion
return 'Diff'(operand, op(var_list));
end proc;