This commit was generated by cvs2svn to compensate for changes in r2,

which included commits to RCS files with non-trunk default branches.


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

5 files changed, 389 insertions, 0 deletions

diff --git a/src/gr/H.c b/src/gr/H.c
new file mode 100644
index 0000000..c345f78
--- /dev/null
+++ b/src/gr/H.c
@@ -0,0 +1,159 @@
+      t1 = g_uu[2][2];
+      t2 = 1/r;
+      t4 = X_ud[0][2];
+      t5 = Diff(h,rho);
+      t7 = X_ud[1][2];
+      t8 = Diff(h,sigma);
+      t10 = zz*t2-t4*t5-t7*t8;
+      t11 = t1*t10;
+      t12 = g_uu[1][1];
+      t14 = X_ud[0][1];
+      t16 = X_ud[1][1];
+      t18 = yy*t2-t14*t5-t16*t8;
+      t19 = t12*t18;
+      t21 = r*r;
+      t23 = 1/r/t21;
+      t30 = Diff(h,rho,rho);
+      t33 = Diff(h,rho,sigma);
+      t38 = Diff(h,sigma,sigma);
+      t40 = -yy*zz*t23-X_udd[0][1][2]*t5-X_udd[1][1][2]*t8-t14*t4*t30-t16*t4*
+t33-t14*t7*t33-t16*t7*t38;
+      t44 = g_uu[0][2];
+      t45 = t44*t44;
+      t46 = t10*t10;
+      t48 = yy*yy;
+      t49 = zz*zz;
+      t56 = X_ud[0][0];
+      t57 = t56*t56;
+      t59 = X_ud[1][0];
+      t63 = t59*t59;
+      t65 = (t48+t49)*t23-X_udd[0][0][0]*t5-X_udd[1][0][0]*t8-t57*t30-2.0*t59*
+t56*t33-t63*t38;
+      t70 = xx*t2-t56*t5-t59*t8;
+      t71 = t70*t70;
+      t72 = t70*t71;
+      t74 = Diff(g_uu[0][0],zz);
+      t77 = g_uu[1][2];
+      t78 = t77*t10;
+      t79 = xx*xx;
+      t86 = t14*t14;
+      t91 = t16*t16;
+      t93 = (t79+t49)*t23-X_udd[0][1][1]*t5-X_udd[1][1][1]*t8-t86*t30-2.0*t16*
+t14*t33-t91*t38;
+      t97 = g_uu[0][1];
+      t99 = Diff(g_uu[0][0],yy);
+      t102 = t77*t18;
+      t109 = t4*t4;
+      t114 = t7*t7;
+      t116 = (t79+t48)*t23-X_udd[0][2][2]*t5-X_udd[1][2][2]*t8-t109*t30-2.0*t7*
+t4*t33-t114*t38;
+      t120 = t44*t10;
+      t135 = -xx*zz*t23-X_udd[0][0][2]*t5-X_udd[1][0][2]*t8-t56*t4*t30-t59*t4*
+t33-t56*t7*t33-t59*t7*t38;
+      t139 = g_uu[0][0];
+      t140 = t139*t70;
+      t141 = Diff(g_uu[1][2],xx);
+      t142 = t141*t10;
+      t145 = t44*t70;
+      t146 = Diff(g_uu[1][2],zz);
+      t147 = t146*t10;
+      t151 = Diff(g_uu[0][0],xx);
+      t154 = t10*t46;
+      t156 = Diff(g_uu[2][2],zz);
+      t159 = t97*t70;
+      t160 = Diff(g_uu[1][2],yy);
+      t161 = t160*t10;
+      t164 = -2.0*t11*t19*t40-t45*t46*t65-t44*t72*t74/2.0-2.0*t78*t19*t93-t97*
+t72*t99/2.0-2.0*t11*t102*t116-2.0*t120*t102*t135-t140*t142*t18-t145*t147*t18-
+t139*t72*t151/2.0-t1*t154*t156/2.0-t159*t161*t18;
+      t166 = Diff(g_uu[2][2],yy);
+      t170 = Diff(g_uu[2][2],xx);
+      t173 = t77*t46;
+      t188 = -xx*yy*t23-X_udd[0][0][1]*t5-X_udd[1][0][1]*t8-t56*t14*t30-t59*t14
+*t33-t56*t16*t33-t59*t16*t38;
+      t192 = t18*t18;
+      t193 = t18*t192;
+      t195 = Diff(g_uu[1][1],zz);
+      t199 = Diff(g_uu[1][1],xx);
+      t202 = t195*t192;
+      t205 = t199*t192;
+      t209 = Diff(g_uu[1][1],yy);
+      t212 = t1*t1;
+      t215 = t1*t46;
+      t216 = t77*t40;
+      t219 = t97*t192;
+      t221 = t12*t192;
+      t223 = -t77*t154*t166/2.0-t44*t154*t170/2.0-2.0*t173*t44*t188-t77*t193*
+t195/2.0-t97*t193*t199/2.0-t145*t202/2.0-t140*t205/2.0-t12*t193*t209/2.0-t212*
+t46*t116-2.0*t215*t216-t219*t142-t221*t161;
+      t225 = t74*t71;
+      t228 = t99*t71;
+      t235 = t151*t71;
+      t238 = t156*t46;
+      t241 = t77*t192;
+      t243 = t209*t192;
+      t246 = t77*t77;
+      t249 = t97*t18;
+      t252 = t170*t46;
+      t255 = t166*t46;
+      t258 = -t102*t225/2.0-t19*t228/2.0-t11*t225/2.0-t78*t228/2.0-t120*t235/
+2.0-t145*t238/2.0-t241*t147-t159*t243/2.0-t246*t46*t93-t249*t235/2.0-t140*t252/
+2.0-t159*t255/2.0;
+      t263 = Diff(g_uu[0][1],zz);
+      t266 = Diff(g_uu[0][1],yy);
+      t267 = t266*t70;
+      t269 = t139*t71;
+      t270 = Diff(g_uu[0][1],xx);
+      t271 = t270*t18;
+      t273 = t97*t71;
+      t274 = Diff(g_uu[0][2],yy);
+      t275 = t274*t10;
+      t277 = Diff(g_uu[0][2],xx);
+      t278 = t277*t10;
+      t282 = t266*t18;
+      t285 = t263*t18;
+      t290 = t44*t71;
+      t291 = Diff(g_uu[0][2],zz);
+      t292 = t291*t10;
+      t297 = -t249*t252/2.0-t102*t238/2.0-t241*t263*t70-t221*t267-t269*t271-
+t273*t275-t269*t278-t19*t255/2.0-t78*t282*t70-t11*t285*t70-t120*t271*t70-t290*
+t292-2.0*t273*t139*t188;
+      t303 = t97*t97;
+      t313 = t97*t188;
+      t316 = t12*t12;
+      t332 = -2.0*t290*t97*t40-t303*t71*t93-2.0*t290*t139*t135-t45*t71*t116-
+t303*t192*t65-2.0*t221*t313-t316*t192*t93-2.0*t241*t97*t135-2.0*t241*t12*t40
+-2.0*t45*t10*t70*t135-t246*t192*t116-t273*t282;
+      t333 = t139*t139;
+      t336 = t145*t40;
+      t351 = t44*t46;
+      t357 = t291*t70;
+      t362 = t159*t40;
+      t365 = -t333*t71*t65-2.0*t78*t336-t249*t278*t70-t102*t292*t70-t19*t275*
+t70-2.0*t11*t249*t135-2.0*t78*t249*t188-t351*t277*t70-t173*t274*t70-t290*t285-
+t215*t357-2.0*t120*t19*t188-2.0*t11*t362;
+      t371 = t140*t65;
+      t374 = t140*t135;
+      t379 = t140*t188;
+      t382 = t159*t93;
+      t401 = -2.0*t246*t10*t18*t40-2.0*t120*t371-2.0*t11*t374-t219*t270*t70-2.0
+*t78*t379-2.0*t78*t382-2.0*t120*t159*t188-2.0*t19*t382-2.0*t303*t18*t70*t188
+-2.0*t19*t379-2.0*t249*t371-2.0*t249*t145*t135;
+      t406 = t145*t116;
+      t422 = t44*t135;
+      t427 = t146*t18;
+      t431 = -2.0*t102*t362-2.0*t102*t374-2.0*t11*t406-2.0*t102*t406-2.0*t19*
+t336-t11*t202/2.0-t78*t243/2.0-2.0*t120*t249*t65-t120*t205/2.0-2.0*t215*t422-
+t351*t141*t18-t215*t427-t173*t160*t18;
+      t441 = t269+2.0*t249*t70+2.0*t120*t70+t221+2.0*t78*t18+t215;
+      t442 = sqrt(t441);
+      t452 = t151*t70+t267+t357+t271+t209*t18+t427+t278+t161+t156*t10+t139*t65+
+2.0*t313+2.0*t422;
+      t456 = partial_d_ln_sqrt_g[0];
+      t459 = partial_d_ln_sqrt_g[1];
+      t462 = partial_d_ln_sqrt_g[2];
+      t477 = t12*t93+2.0*t216+t1*t116+t139*t456*t70+t97*t459*t70+t44*t462*t70+
+t97*t456*t18+t12*t459*t18+t77*t462*t18+t44*t456*t10+t77*t459*t10+t1*t462*t10;
+      H = (t164+t223+t258+t297+t332+t365+t401+t431)/t442/t441+(t452+t477)/t442+
+(K_uu[0][0]*t71+2.0*K_uu[0][1]*t18*t70+2.0*K_uu[0][2]*t10*t70+K_uu[1][1]*t192+
+2.0*K_uu[1][2]*t10*t18+K_uu[2][2]*t46)/t441-K;
diff --git a/src/gr/Makefile b/src/gr/Makefile
new file mode 100644
index 0000000..892450c
--- /dev/null
+++ b/src/gr/Makefile
@@ -0,0 +1,29 @@
+# Makefile for GR code
+# $Id: Makefile,v 1.1.1.1 2001-06-14 11:53:21 jthorn Exp $
+
+#
+# Environment Variables:
+#	MAPLE_VERSION	used via @ifdef for version control in Maple code;
+#			typically set to something like  MAPLE_V_RELEASE_4
+#
+# Targets:
+#	mm		==> preprocess all *.maple files to produce *.mm files
+#	clean		==> delete *.mm (mpp Maple preprocessor output)
+#	run		==> load all the code into Maple
+#
+
+###############################################################################
+
+.PHONY	: mm
+mm	: $(patsubst %.maple, %.mm, $(wildcard *.maple))
+
+%.mm	: %.maple $(wildcard *.minc $(gfa_dir)/*.minc $(gfa_dir)/*.maple)
+	mpp -D$(MAPLE_VERSION) <$< >$@
+
+.PHONY	: clean
+clean	:
+	-rm *.mm
+
+.PHONY	: run
+run	: mm
+	maple -f <setup.mm >maple.log
diff --git a/src/gr/gr_gfas.minc b/src/gr/gr_gfas.minc
new file mode 100644
index 0000000..92e781e
--- /dev/null
+++ b/src/gr/gr_gfas.minc
@@ -0,0 +1,25 @@
+# Maple include file listing all gfa functional-dependence forms
+
+g_dd,
+K_dd,
+
+g_uu, g_uu__fnd,
+K_uu, K_uu__fnd,
+
+K, K__fnd,
+
+partial_d_ln_sqrt_g, partial_d_ln_sqrt_g__fnd,
+
+h, h__fnd,
+
+s_d, s_d__fnd,
+partial_d_s_d, partial_d_s_d__fnd,
+
+n_u, n_u__fnd,
+
+H, H__fnd,
+
+HA, HA__fnd,
+HB, HB__fnd,
+HC, HC__fnd,
+HD, HD__fnd # no comma
diff --git a/src/gr/horizon_function.maple b/src/gr/horizon_function.maple
new file mode 100644
index 0000000..c1750ea
--- /dev/null
+++ b/src/gr/horizon_function.maple
@@ -0,0 +1,127 @@
+# horizon.maple - evaluate LHS of horizon function
+#
+# non_unit_normal - compute s_d
+# non_unit_normal_gradient - compute partial_d_s_d
+# horizon_function - compute HA,HB,HC,HD,H = fn(s_d)
+#
+
+################################################################################
+################################################################################
+################################################################################
+
+#
+# This function computes the non-unit outward-pointing normal vector
+# (field) to the horizon, using the equation on p.7.1 of my apparent
+# horizon finding notes.
+#
+non_unit_normal :=
+proc()
+global
+       @include "../maple/coords.minc",
+       @include "../gr/gr_gfas.minc";
+local i, u;
+
+assert_fnd_exists(h);
+assert_fnd_exists(X_ud);
+assert_fnd_exists(s_d, fnd);
+
+	for i from 1 to N
+	do
+	s_d__fnd[i] := x_xyz[i]/r
+		       - sum('X_ud[u,i]*Diff(h,x_rrs[u])', 'u'=A_min..N);
+	end do;
+
+NULL;
+end proc;
+
+
+################################################################################
+
+#
+# This function computes the gradient of the non-unit outward-pointing
+# normal vector (field) to the horizon, using the equation on p.7.2 of
+# my apparent horizon finding notes.
+#
+non_unit_normal_gradient :=
+proc()
+global
+       @include "../maple/coords.minc",
+       @include "../gr/gr_gfas.minc";
+local temp, i, j, k, u, v;
+
+assert_fnd_exists(h);
+assert_fnd_exists(X_ud);
+assert_fnd_exists(X_udd);
+assert_fnd_exists(partial_d_s_d, fnd);
+
+	for i from 1 to N
+	do
+	for j from 1 to N
+	do
+	temp := `if`(i=j, sum('x_xyz[k]^2','k'=1..3), 0);
+	# simplify likes to put things over a common denominator,
+	# so we only apply it to the individual terms, not to
+	# the whole experssion
+	partial_d_s_d__fnd[i,j]
+	   := + simplify(  (temp-x_xyz[i]*x_xyz[j])/r^3  )
+	      - simplify(  sum('X_udd[u,i,j]*Diff(h,x_rrs[u])', 'u'=A_min..N)  )
+	      - simplify(  msum('X_ud[u,i]*X_ud[v,j]*Diff(h,x_rrs[u],x_rrs[v])',
+				'u'=A_min..N, 'v'=A_min..N)  );
+	end do;
+	end do;
+
+NULL;
+end proc;
+
+################################################################################
+
+#
+# This function computes the LHS of the apparent horizon equation
+# as a function of 1st and 2nd angular partial derivatives of the
+# apparent horizon radius $h$, in the Maple gridfn array H.fnd ,
+# using equation (15) in my 1996 apparent horizon finding paper.
+#
+horizon_function :=
+proc()
+global
+       @include "../maple/coords.minc",
+       @include "../gr/gr_gfas.minc";
+local i,j,k,l;
+
+assert_fnd_exists(g_uu);
+assert_fnd_exists(K_uu);
+assert_fnd_exists(partial_d_ln_sqrt_g);
+assert_fnd_exists(s_d, fnd);
+assert_fnd_exists(partial_d_s_d, fnd);
+
+#
+# n.b. no simplify() in these subexpressions; empirically it makes
+#      them much *more* complicated (measured by codegen(...,optimized))
+#
+
+# (15a) in my paper
+HA__fnd := - msum('g_uu[i,k]*s_d__fnd[k] * g_uu[j,l]*s_d__fnd[l]
+					 * partial_d_s_d__fnd[i,j]',
+		  'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N)
+	   - (1/2)*msum('g_uu[i,j]*s_d__fnd[j]
+			* Diff(g_uu[k,l],x_xyz[i]) * s_d__fnd[k]*s_d__fnd[l]',
+			'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N);
+
+# (15b) in my paper
+HB__fnd := + msum('Diff(g_uu[i,j],x_xyz[i])*s_d__fnd[j]', 'i'=1..N, 'j'=1..N)
+	   + msum('g_uu[i,j]*partial_d_s_d__fnd[i,j]', 'i'=1..N, 'j'=1..N)
+	   + msum('g_uu[i,j]*partial_d_ln_sqrt_g[i]*s_d__fnd[j]',
+		  'i'=1..N, 'j'=1..N);
+
+# (15c) in my paper
+HC__fnd := msum('K_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N);
+
+# (15d) in my paper
+HD__fnd := msum('g_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N);
+
+# n.b. no simplify() here, it would try to put things over a common
+#      denominator, which would make the equation much more complicated
+H__fnd := HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K;
+
+NULL;
+end proc;
diff --git a/src/gr/setup_gfas.maple b/src/gr/setup_gfas.maple
new file mode 100644
index 0000000..5da990b
--- /dev/null
+++ b/src/gr/setup_gfas.maple
@@ -0,0 +1,49 @@
+# Maple code to set up all gridfn arrays
+# $Id$
+
+#
+# setup_gr_gfas - setup all the GR gfas
+#
+
+################################################################################
+
+setup_gr_gfas :=
+proc()
+global
+       @include "../gr/gr_gfas.minc";
+
+# basic metric & extrinsic curvature fields
+make_gfa('g_dd', {inert}, [1..N, 1..N], symmetric);
+make_gfa('K_dd', {inert}, [1..N, 1..N], symmetric);
+
+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
+make_gfa('partial_d_ln_sqrt_g', {inert,fnd}, [1..N], none);
+
+# radius of horizon
+make_gfa('h', {inert, fnd}, [], none);
+
+# outward-pointing *non*-unit normal 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);
+
+# outward-pointing unit normal to horizon
+make_gfa('n_u', {inert,fnd}, [1..N], none);
+
+# LHS of apparent horizon equation
+make_gfa('H', {inert, fnd}, [], 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);
+
+NULL;
+end proc;