diff options
author | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2001-06-14 11:53:21 +0000 |
---|---|---|
committer | jthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5> | 2001-06-14 11:53:21 +0000 |
commit | 6a0395d7e63372be04d20ee447b5d88af4f09f69 (patch) | |
tree | 925ce0fd4046ff4bf7b05f7f1c7672be858836fe | |
parent | 54ea67b814e9e5e8a00c797b320c0ad2b02bf4f0 (diff) |
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
-rw-r--r-- | README | 51 | ||||
-rw-r--r-- | interface.ccl | 2 | ||||
-rw-r--r-- | param.ccl | 2 | ||||
-rw-r--r-- | schedule.ccl | 2 | ||||
-rw-r--r-- | src/gr/H.c | 159 | ||||
-rw-r--r-- | src/gr/Makefile | 29 | ||||
-rw-r--r-- | src/gr/gr_gfas.minc | 25 | ||||
-rw-r--r-- | src/gr/horizon_function.maple | 127 | ||||
-rw-r--r-- | src/gr/setup_gfas.maple | 49 | ||||
-rw-r--r-- | src/make.code.defn | 9 | ||||
-rw-r--r-- | src/maple/Diff.maple | 260 | ||||
-rw-r--r-- | src/maple/Makefile | 29 | ||||
-rw-r--r-- | src/maple/coeffs.maple | 10 | ||||
-rw-r--r-- | src/maple/coords.maple | 68 | ||||
-rw-r--r-- | src/maple/coords.minc | 18 | ||||
-rw-r--r-- | src/maple/gfa.maple | 315 | ||||
-rw-r--r-- | src/maple/setup.maple | 24 | ||||
-rw-r--r-- | src/maple/util.maple | 234 |
18 files changed, 1413 insertions, 0 deletions
@@ -0,0 +1,51 @@ +Cactus Code Thorn AHFinderDirect +Authors : Jonathan Thornburg <jthorn@aei.mpg.de> +CVS info : $Header$ +-------------------------------------------------------------------------- + +Purpose of the thorn +==================== + +This thorn finds an apparent horizon given the 3D xyz-grid metric and +extrinsic curvature. It uses a direct method, writing the apparent +horizon equation as an elliptic PDE on angular-coordinate space. This +is fast, but does require an intitial guess for the apparent horizon +position. + + +Documentation +============= + +The doc/ directory contains a user's guide and a programmer's guide +for this thorn. + + +Distribution +============ + +This thorn should eventually become part of the CactusEinstein arrangement. + +This thorn is distributed under terms of the GNU public license; see +the file COPYING in this directory for details. + + +Other Software Required +======================= + +* This thorn requires the Einstein thorn. +* This thorn is written in C++, so you'll need a C++ compiler to compile + this thorn. It calls various C and Fortran numerical libraries; if + these aren't already installed on your system, then you'll need C + and Fortran compilers to compile them. +* Most of this thorn's relativity code is machine-generated using Maple + (version 6), but you don't need Maple unless you want to modify the + relativity code. + + +Portability +=========== + +This thorn is known to compile and run successfully on the following +systems: +* Red Hat GNU/Linux, gcc 2.95.3 +* Alpha OSF1 V4.0, gcc 2.95.3 diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..3905af4 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,2 @@ +# Interface definition for thorn AHFinderDirect +# $Header$ diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..9e673cc --- /dev/null +++ b/param.ccl @@ -0,0 +1,2 @@ +# Parameter definitions for thorn AHFinderDirect +# $Header$ diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..815801e --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,2 @@ +# Schedule definitions for thorn AHFinderDirect +# $Header$ 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; diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..a3d158e --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,9 @@ +# Main make.code.defn file for thorn AHFinderDirect +# $Header$ + +# Source files in this directory +SRCS = + +# Subdirectories containing source files +SUBDIRS = + diff --git a/src/maple/Diff.maple b/src/maple/Diff.maple new file mode 100644 index 0000000..c21d79b --- /dev/null +++ b/src/maple/Diff.maple @@ -0,0 +1,260 @@ +# Diff.maple -- Diff() and friends +# $Id$ + +# +# simplify/Diff - simplify(...) function for expressions containing Diff() +# Diff - "inert" derivative function (low-level simplifications) +# Diff/gridfn - simplifications that know about gridfns +# + +############################################################################### +############################################################################### +############################################################################### + +# +# This function is called automatically by simplify() whenever +# simplify()'s argument contains a Diff() call. This interface +# to simplify() is described in +# Michael B. Monagan +# "Tips for Maple Users" +# The Maple Technical Newsletter +# Issue 10 (Fall 1993), p.10-18, especially p.17-18 +# but not in the Maple library manual or the online help, which both +# imply that simplify(..., Diff) is needed for this procedure to be +# invoked). +# +# This function recurses over the expression structure. For each +# Diff() call, it calls Diff() itself to do "basic" derivative +# simplifications, then calls `Diff/gridfn`() to do further simplifications +# based on gridfn properties. +# +# expr = (in) An expression to be simplified. +# +# Results: +# The function returns the "simplified" expression. +# +`simplify/Diff` := +proc(expr) +option remember; # performance optimization +local temp; + +# recurse on tables, equations, lists, sets, sums, products, powers +if type(expr, {table, `=`, list, set, `+`, `*`, `^`}) + then return map(simplify, expr); +end if; + +# return unchanged on numbers, names, procedures +if type(expr, {numeric, name, procedure}) + then return expr; +end if; + +# function call ==> simplify Diff() calls, recurse on others +if type(expr, function) + then + if (op(0, expr) = 'Diff') + then return Diff(op(expr)); # "basic" derivative simplifications + else + temp := map(simplify, [op(expr)]); + op(0,expr); return '%'(op(temp)); + end if; +end if; + +# unknown type +ERROR(`expr has unknown type`, `whattype(expr)=`, whattype(expr)); + +end proc; + +############################################################################### + +# +# This function implements some of the more useful simplification +# rules of algebra and calculus for Maple's "inert" differentiation +# function Diff() . It then calls `Diff/gridfn`() to do any +# further simplifications based on gridfn properties. +# +# In particular: +# - It distributes over sums. +# - It distributes over leading numeric factors. (Maple's built-in +# simplifier will always make numeric factors leading.) +# - It knows that derivatives of numeric constants are zero. +# - It knows that the derivative of a name with respect to itself is one. +# - It knows the product rule +# Diff(f*g,x) = Diff(f,x)*g + f*Diff(g,x) +# - It knows the power rule +# Diff(f^n,x) = n * f^(n-1) * Diff(f,x) +# - It flattens multi-level "Diff() trees", i.e. it makes the transformation +# Diff(Diff(f,x),y) = Diff(f,x,y); +# - It canonicalizes the order of the things we're differentiating with +# respect to, so that simplify() will recognize that derivatives +# commute, i.e. that +# Diff(f,x,y) = Diff(f,y,x) +# +# Arguments: +# operand = (in) The thing to be differentiated. +# var_seq = (in) (varargs) An expression sequence of the variables to +# differentiate with respect to. +# +unprotect('Diff'); +Diff := +proc(operand) # varargs +option remember; # performance optimization +local var_list, + nn, nderiv, + op_cdr, + abs_operand, + f, g, x, x_car, x_cdr, temp, + n, + inner_operand, inner_var_list, k, + sorted_var_list; + +var_list := [args[2..nargs]]; + +nn := nops(operand); +nderiv := nops(var_list); + +##print(`Diff(): nn=`,nn,` nderiv=`,nderiv, +## ` operand=`,operand,` var_list=`,var_list); + +# distribute (recurse) over sums +if (type(operand, `+`)) + then return simplify( + sum('Diff(op(k, operand), op(var_list))', 'k'=1..nn) + ); +end if; + +# distribute (recurse) over leading numeric factors +if (type(operand, `*`) and (nn >= 1) and type(op(1,operand), numeric)) + then + op_cdr := product('op(k,operand)', 'k'=2..nn); + return simplify( op(1,operand) * Diff(op_cdr, op(var_list)) ); +end if; + +# derivatives of numbers are zero +if (type(operand, numeric)) + then return 0; +end if; + +# derivative of a name with respect to itself is one +if (type(operand, name) and (var_list = [operand])) + then return 1; +end if; + +# implement the product rule +# ... we actually implement it only for a 1st derivative of the +# product of two factors; for fancier cases we recurse +if (type(operand, `*`)) + then + if (nn = 0) + then # empty product = 1 + return 0; # ==> Diff(1) = 0 + + elif (nn = 1) + then # singleton product is equal to the single factor + return simplify( Diff(op(1,operand), op(var_list)) ); + + elif (nn >= 2) + then # nontrivial product rule case + + # note that if we have more than 2 factors, g will + # be the product of the cdr of the factor list, so + # our later Diff(g,...) calls will recurse as + # necessary to handle this case + + f := op(1,operand); + g := product('op(k,operand)', 'k'=2..nn); + + if (nderiv = 1) + then # product rule for 1st derivatives: + # Diff(f*g,x) = Diff(f,x)*g + f*Diff(g,x) + x := var_list[1]; + + return simplify( Diff(f,x)*g + f*Diff(g,x) ); + + else # product rule for 2nd (and higher) derivatives + # ==> recurse on derivatives + x_car := var_list[1]; + x_cdr := var_list[2..nderiv]; + temp := simplify( Diff(f*g, x_car) ); + return simplify( Diff(temp, op(x_cdr)) ); + end if; + else ERROR(`impossible value of nn in product rule!`, + `(this should never happen!)`, + ` operand=`,operand,` nn=`,nn); + end if; +end if; + +# implement the power rule +# ... we actually implement it only for 1st derivatives; for higher +# derivatives we recurse +if (type(operand, `^`)) + then + f := op(1,operand); + n := op(2,operand); + + if (nderiv = 1) + then # power rule for 1st derivatives: + # Diff(f^n,x) = n * f^(n-1) * Diff(f,x) + x := var_list[1]; + + return simplify( n * f^(n-1) * Diff(f, x) ); + + else # power rule for 2nd (and higher) derivatives: + # ==> recurse on derivatives + x_car := var_list[1]; + x_cdr := var_list[2..nderiv]; + temp := simplify( Diff(f^n, x_car) ); + + return simplify( Diff(temp, op(x_cdr)) ); + end if; +end if; + +# flatten (recurse on) multi-level "Diff() trees" +if (type(operand, function) and (op(0, operand) = 'Diff')) + then + inner_operand := op(1, operand); + inner_var_list := [op(2..nn, operand)]; + + return simplify( + Diff(inner_operand, op(inner_var_list), op(var_list)) + ); +end if; + +# get to here +# ==> no other transformations to make +# ==> canonicalize ordering of variables and return reconstructed Diff() call +sorted_var_list := sort_var_list(var_list); +return `Diff/gridfn`(operand, op(sorted_var_list)); + +end proc; + +############################################################################### + +# +# This function implements further simplification rules for Diff() +# based on gridfn properties. +# +# It currently knows about the following simplifications: +# - Diff(X_ud[u,i], x_xyz[j]) --> X_udd[u,i,j] +# +# Arguments: +# operand = (in) The thing to be differentiated. +# var_seq = (in) (varargs) An expression sequence of the variables to +# differentiate with respect to. +# +`Diff/gridfn` := +proc(operand) # varargs +option remember; # performance optimization +global + @include "../maple/coords.minc"; +local var_list, posn; + +var_list := [args[2..nargs]]; + +if ( type(operand, indexed) and (op(0,operand) = 'X_ud') + and (nops(var_list) = 1) and member(var_list[1],x_xyz_list,'posn') ) + then return X_udd[op(operand), posn]; + else return 'Diff'(operand, op(var_list)); # return is *unevaluated* + # to avoid infinite recursion! +end if; + +end proc; diff --git a/src/maple/Makefile b/src/maple/Makefile new file mode 100644 index 0000000..c46eea4 --- /dev/null +++ b/src/maple/Makefile @@ -0,0 +1,29 @@ +# Makefile for maple 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/maple/coeffs.maple b/src/maple/coeffs.maple new file mode 100644 index 0000000..b45114b --- /dev/null +++ b/src/maple/coeffs.maple @@ -0,0 +1,10 @@ +# coeffs.maple - set up coordinate coefficients +# $Id$ +# +# setup_X_ud +# setup_X_udd +# + +################################################################################ + + diff --git a/src/maple/coords.maple b/src/maple/coords.maple new file mode 100644 index 0000000..4e8953d --- /dev/null +++ b/src/maple/coords.maple @@ -0,0 +1,68 @@ +# coords.maple -- set up and covert between coordinates, gridfn arrays, etc +# $Id$ + +# +# setup_coords - set up coordinates +# setup_coord_gfas - setup up coordinate-transformation gfas +# + +############################################################################### + +# +# This function sets up the coordinates. +# It also sets up the Kronecker delta +# +setup_coords := +proc() +global +@include "coords.minc"; + +N := 3; # number of spatial coordinates, + # i.e. spatial coordinate indices + # are 1..N +N_ang := 2; # number of angular spatial coordinates +A_min := 2; # min angular-spatial-coordinate + # spatial coordinate index, + # i.e. angular spatial coordinate + # indices are A_min..N + +delta := array(1..N, 1..N, identity); # Kronecker delta + +# (x,y,z) coordinates +x_xyz := array(1..N); +x_xyz[1] := xx; +x_xyz[2] := yy; +x_xyz[3] := zz; +x_xyz_list := map(op, [entries(x_xyz)]); + +# (r,rho,sigma) coordinates +# we hard-wire that r is the first coordinate and remaining ones are angular +x_rrs := array(1..N); +x_rrs[1] := r; +x_rrs[2] := rho; +x_rrs[3] := sigma; +x_rrs_list := map(op, [entries(x_rrs)]); + +# list of all coordinates (for sorting things into "nice" order) +x_all_list := [op(x_xyz_list), op(x_rrs_list)]; + +NULL; +end proc; + +############################################################################### + +# +# This function sets up the coordinate-transformation gfas +# +setup_coord_gfas := +proc() +global +@include "coords.minc"; + +# 1st ($X^u_i$) and 2nd ($X^u_{ij}$) partial derivs +# of (r,rho,sigma) wrt (xx,yy,zz) +make_gfa('X_ud', {inert,fnd}, [A_min..N, 1..N], none); +make_gfa('X_udd', {inert,fnd}, [A_min..N, 1..N, 1..N], symmetric3_23); + +NULL; +end proc; diff --git a/src/maple/coords.minc b/src/maple/coords.minc new file mode 100644 index 0000000..532bfa0 --- /dev/null +++ b/src/maple/coords.minc @@ -0,0 +1,18 @@ +# coords.minc -- global variables for coordinates +# $Id$ + +# coordinates +N, # number of spatial coordinates +N_ang, # number of angular spatial coordinates +A_min, # min angular-spatial-coordinate coordinate index +delta, # Kronecker delta +x_xyz, x_xyz_list, # array and list of (xx,yy,zz) coordinates +x_rrs, x_rrs_list, # array and list of (r,rho,sigma) coordinates +r, rho, sigma, # individual patch coordinates +xx, yy, zz, # Cartesian coordinates +x_all_list, # list of all coordinates + +# coordinate-transformation gfas +X_ud, X_ud__fnd, +X_udd, X_udd_fnd + # no comma diff --git a/src/maple/gfa.maple b/src/maple/gfa.maple new file mode 100644 index 0000000..e0733b6 --- /dev/null +++ b/src/maple/gfa.maple @@ -0,0 +1,315 @@ +# gfa.maple -- low-level gridfn array "database" +# $Id$ +# +# ***** gridfn arrays ***** +# ***** functional dependences ***** +# ***** gridfn-array sets ***** +# make_gfa - construct a new gridfn array +# assert_fnd_exists - assert (verify) that a given gridfn array fnd form exists +# fnd_name - Maple name of functional-dependence form +# +# index/symmetric3_23 - Maple indexing fn for rank 3 array sym on axes 23 +# +# print_symmetric3_23 - prettyprint a symmetric3_23 array +# + +############################################################################### + +# +# ***** gridfn arrays ***** +# + +# +# The basic data objects for the Maple code are grid function arrays, +# a.k.a. gridfn arrays, such as the metric g_dd, the Christoffel symbols +# Gamma_udd, etc etc. +# +# By convention, the name of a nonscalar geometric object ends with an +# underscore ("_") followed by a sequence of "d" or "u" specifying whether +# the indices are down (covariant) or up (contravariant). Thus, for +# example, g_dd is the covariant metric, while g_uu is the inverse metric. +# +# Scalar geometric objects are represented by Maple scalars, while +# nonscalar objects (often but not necessarily tensors) are represented +# by Maple arrays, using Maple indexing functions for any symmetries +# or sparsity. For example, the metric g_dd is represented by an +# array +# array(1..N, 1..N, symmetric) +# By convention, for a symmetric array of this type, we only compute +# the upper triangle. +# + +############################################################################### + +# +# ***** functional dependences ***** +# + +# +# Within the Maple code, we distinguish between two distinct forms +# that a gridfn array may take: +# - The "inert" form of a gridfn array represents (only) the symmetries, +# zero (or nonzero but constant) elements, and indexing semantics of +# the gridfn array. +# - The "functional dependence" form of a gridfn array additionally +# represents how that gridfn array is computed from other gridfn +# arrays. +# +# For example, consider the inverse metric g_uu . It's inert form +# is a 3*3 symmetric array, with all elements being unassigned symbols +# (i.e. evaluating to themselves, as is usual in Maple). It's functional +# dependence form, in contrast, is also a 3*3 symmetric array, but its +# elements are expressions giving the inverse-metric elements in terms +# of components of the metric and/or the metric determinant. +# +# Any given gridfn array may exist in either or (the usual case) both +# an inert form and a functional dependence forms. There may also be +# multiple functional dependence forms, corresponding to alternative +# formulas for its computation.) +# +# The inert form of a gridfn array is represented by a Maple variable +# or array (see below) of the same name (eg g_uu ), unassigned except +# for any symmetries and/or zero components the gridfn array has. +# +# The functional dependence form of a gridfn array is represented by +# a Maple array (see below) with a name formed by appending "__fnd" +# to the gridfn array's name (eg g_uu__fnd ); this array has the same +# symmetries as the corresponding inert Maple array, but only those +# elements needed will be assigned. If there are multiple functional +# dependence forms, by convention each uses a name beginning in "__" +# and ending in "fnd", eg H__ps_fnd . +# +# The PDE compiler normally generates C code referring to the (inert) +# gridfn array names, so these should be valid C variable names. +# + +# +# To allow some run-time error checking, we keep a global table +# of all the gridfn array functional-dependence forms: +# `gfa/fnd_table` +# ... table index is gfa__fnd_name +# ... table value is true +# assert_fnd_exists() uses this table to verify that a given gridfn +# array functional-dependence form does indeed exist. +# + +# +# Functional-dependence forms may be multilevel, as in (eg) +# g_dd__Kerr_fnd__ps_fnd . cat() may be used to combine names in +# this manner. +# + +############################################################################### + +# +# This function makes the Maple representation of a single gridfn array. +# +# Arguments: +# gfa_name = (out*) The Maple variable to be created for the gridfn +# array. This argument should be passed as an +# unevaluated (quoted) name. +# fnd_name_set = (in) A set of functional dependence types the gridfn +# array is to have. The set's members should be +# drawn from from {inert, fnd, ...} as appropriate, +# where the special value inert refers to the +# inert form. +# index_bounds = (in) The list of index bounds for the gridfn array, +# index_fn = (in) The name of the array indexing function, if any, or +# 'none' if the gridfn array has no symmetries and thus +# should be represented by a Maple array with the standard +# "ordinary array" indexing semantics. +# +# Thus, for example, the calls +# make_gfa('g', {inert, fnd}, [], none); +# make_gfa('g_dd', {inert}, [1..N, 1..N], symmetric); +# make_gfa('R_dd', {inert, fnd}, [1..N, 1..N], symmetric); +# make_gfa('Gamma_udd', {inert, fnd}, [1..N, 1..N, 1..N], symmetric3_23); +# would set up the Maple representations of the metric determinant, +# the metric, the Ricci tensor, and the Christoffel symbols respectively. +# +make_gfa := +proc(gfa_name::name, + fnd_name_set::set(name), + index_bounds::list(integer..integer), + index_fn::{identical(none), name}) # note "name" here, not "procedure", + # since we'll get the name of the + # indexing fn, which may not have + # been assigned its procedure yet +global N, `gfa/fnd_table`; +local fnd_name, var_name, index_fn_seq; + +# firewall checks +if (not type(N, integer)) + then ERROR("must set up coordinates first!"); +end if; +if ((index_bounds = []) and (index_fn <> 'none')) + then ERROR("not meaningful to specify a symmetry function for a scalar!", + " gfa_name=",gfa_name); +end if; + +# create the Maple variable(s) + for fnd_name in fnd_name_set + do + var_name := Maple_name(gfa_name, fnd_name); + if (assigned(`gfa/fnd_table`[eval(var_name,1)])) + then ERROR("duplicate gfa/fnd definition!", + " gfa_name=", gfa_name, + " fnd_name_set=", fnd_name_set); + end if; + `gfa/fnd_table`[eval(var_name,1)] := true; + if (index_bounds = []) + then # scalar + unassign( eval(var_name,1) ); + else # array + if (index_fn = 'none') + then index_fn_seq := NULL; + else index_fn_seq := index_fn; + end if; + assign( + eval(var_name,1) = array(index_fn_seq, op(index_bounds)) + ); + end if; + end do; + +NULL; +end proc; + +############################################################################### + +# +# This function checks that a given gridfn array functional-dependence +# form exists. If so, this function is a no-op; if not, this function +# calls ERROR(). Hence after calling this function, the caller can assume +# that the given functional-dependence form does indeed exist. +# +# This function may be called with either 1 or 2 arguments (the two +# types of calls are equivalent). +# +# Arguments (1-argument form): +# gfa_name = The name of the fnd-form to verify, eg 'g_dd.fnd'. +# If this isn't inert, it should be quoted so this function +# gets the name, not whatever value may have been assigned to it. +# +# Arguments (2-argument form): +# gfa_name = The name of the gridfn array itself, eg `g_dd`. +# fnd_name = The name of the fnd itself, eg `fnd`. +# +assert_fnd_exists := +proc(gfa_name::name, fnd_name::name) +global `gfa/fnd_table`; +local var_name; + +var_name := Maple_name(gfa_name, args[2..nargs]); + +if (not assigned(`gfa/fnd_table`[eval(var_name,1)])) + then ERROR("functional-dependence form doesn't exist!", + "var_name=",var_name); +end if; +end proc; + +############################################################################### + +# +# This function computes the Maple name corresponding to a given +# gridfn array and functional dependence form. +# +# Arguments: +# gfa_name = The name of the gridfn array itself, eg `g_dd`. +# fnd_name = (optional) +# The name of the functional-dependence form, eg `fnd`, +# or the special name 'inert'. If this argument is omitted +# it defaults to 'inert'. +# +Maple_name := +proc(gfa_name::name, fnd_name::name) +if ((nargs = 1) or (fnd_name = 'inert')) + then return gfa_name; + else return cat(gfa_name, "__", fnd_name); +end if; +end proc; + +############################################################################### +############################################################################### +############################################################################### + +# +# This function is a Maple indexing function for a rank 3 table/array +# which is symmetric in the 2nd and 3rd indices. This function is +# called automagically by Maple whenever a component of such an array +# is referenced. +# +# The ordering is done using the Maple sort() function, which sorts +# numbers numerically, strings lexicographically, and "other things" +# by machine address. +# +# Arguments: +# ilist = (in) The indices to be used. There must be exactly 3 of these. +# tab = (in out) The table/array to be indexed. This uses any remaining +# indexing methods. +# vlist = (optional) +# If present, this is a list of the values to be assigned to +# the (lvalue) table entry. +# +`index/symmetric3_23` := +proc(ilist::list, tab::table, vlist::list) +local k, i, j, index_seq; + +if (not (nops(ilist) = 3)) + then ERROR(`must have exactly 3 indices!`); +end if; + +# ... note we need eval() here to get actual integers, +# rather than variable names whose values are integers +k := eval(ilist[1]); +i := eval(ilist[2]); +j := eval(ilist[3]); + +index_seq := k, op(sort([i,j])); + +if (nargs = 2) + then return tab[index_seq]; # rvalue + else tab[index_seq] := op(vlist); # lvalue +end if; +end proc; + +############################################################################### +############################################################################### +############################################################################### + +# +# This function prettyprints a symmetric3_23 array +# +# Arguments: +# A = (in) The array to be prettyprinted +# +# Results: +# There is no explicit result; the array is prettyprinted as a side effect. +# +print_symmetric3_23 := +proc(A::array) +local bounds, k, i, j, M23; + +if (op(1, eval(A)) <> 'symmetric3_23') + then ERROR(`can only print symmetric3_23 arrays`); +end if; + +bounds := op(2, eval(A)); + +M23 := array(bounds[2..3], symmetric); + + for k in $(bounds[1]) do + + for i in $(bounds[2]) do + for j in $(bounds[3]) do + if (i <= j) # only need upper triangle in (i,j) + then M23[i,j] := A[k,i,j]; + end if; + end do; + end do; + + printf("[%d] = \n", k); + print(M23); + end do; + +NULL; +end proc; diff --git a/src/maple/setup.maple b/src/maple/setup.maple new file mode 100644 index 0000000..8602d4e --- /dev/null +++ b/src/maple/setup.maple @@ -0,0 +1,24 @@ +# setup.maple -- top-level Maple input file to set up PDE compiler system +# $Id$ + +# +# Note not all files here actually *use* the preprocessor, but for +# convenience we preprocess everything (this way we don't have to keep +# track of which files should be used as .maple and which as .mm). +# + +# +# Note paths are evaluated relative to where read command gets executed, +# not relative to this directory, so we use ../maple/ to try and avoid +# ambiguities... +# + +read "../maple/util.mm"; + +read "../maple/coords.mm"; +setup_coords(); + +read "../maple/Diff.mm"; + +read "../maple/gfa.mm"; +setup_coord_gfas(); diff --git a/src/maple/util.maple b/src/maple/util.maple new file mode 100644 index 0000000..2e574df --- /dev/null +++ b/src/maple/util.maple @@ -0,0 +1,234 @@ +# util.maple -- Maple utility functions +# $Id$ + +# +# msum - multiple-index sum() +# arctan_xy - 2-argument arctan() function +# +# sort_var_list - sort a list of variables into a canonical order +# lexorder_vars - lexorder() for "interesting" variables +# +# gensym - generate a "new" symbol name +# gensym/init - reinitialize the gensym() global naming counter +# +# saveit - optionally save variables for debugging +# + +############################################################################### +############################################################################### +############################################################################### + +# +# This function is identical to the standard Maple library function sum() +# except that it allows multiple summation indices. It's usage is best +# illustrated by a series of examples: +# +# msum('f', 'k'=m..n) +# = sum('f', 'k'=m..n) +# +# msum('f', 'k1'=m1..n1, 'k2'=m2..n2) +# = sum( 'sum('f', 'k1'=m1..n1)' , 'k2'=m2..n2 ) +# +# msum('f', 'k1'=m1..n1, 'k2'=m2..n2, 'k3'=m3..n3) +# = sum( +# 'sum( 'sum('f', 'k1'=m1..n1)' , 'k2'=m2..n2 )' +# , +# 'k3'=m3..n3 +# ) +# +# Note, in particular, that the summation indices are nested with the +# innermost summation appearing first. This is the same nesting convention +# as Fortran I/O statement's implicit do loops use. +# +# Bugs: +# - The quoting in the code is pretty obscure... +# +msum := +proc(fn::algebraic) # varargs +local expr, sum_index; + +if (nargs < 2) + then ERROR("must have two or more arguments") +end if; + +expr := fn; + + for sum_index in [args[2..nargs]] + do + # loop invariant: + # msum(fn, <args processed so far>) = sum(expr, sum_index) + expr; sum_index; expr:= 'sum'(''%%'', ''%''); + end do; + +return eval(expr); +end proc; + +############################################################################### + +# +# This function computes the 4-quadrant arctangent of arctan(y/x) . +# It's "just" a wrapper around the Maple 2-argument arctan function, +# providing the what-I-find-natural (x,y) argument +# +# Arguments: +# (x,y) = (In) The x,y coordinates. +# +arctan_xy := +proc(x::algebraic, y::algebraic) +arctan(y,x); # !!! reversedargument order !!! +end proc; + +################################################################################ +################################################################################ +################################################################################ + +# +# This function sorts a list of names (in practice, a list of coordinates) +# into a canonical order. +# +sort_var_list := +proc(var_list::list(name)) +option remember; # performance optimization + +# only get to here the first time we see a given variable list +# (i.e. if it's not in the remember table) +RETURN( sort(var_list, lexorder_vars) ); +end; + +################################################################################ + +# +# This function defines a lexical ordering on variable names. +# +# Arguments: +# x, y = Two names to be compared. +# +# Results: +# This function returns true iff a < b, false otherwise. +# +lexorder_vars := +proc(x::name, y::name) +option remember; # performance optimization +global + @include "coords.minc"; + +if (member(x, x_all_list, 'xposn') and member(y, x_all_list, 'yposn')) + then return evalb(xposn < yposn); + else return lexorder(x, y); +fi; + +end; + +############################################################################### +############################################################################### +############################################################################### + +# +# This procedure generates a (presumably) new symbol name of the +# form +# cat(base_name,count) +# where base_name is passed as an (optional, with a default value +# supplied if it's omitted) argument, and count is the current +# value of the global variable +# `gensym/counter` +# which is postincremented after being used. If the counter is +# unassigned on entry to this function, `gensym/init`() is first +# called to initialize it. +# +# `gensym/init`() may be used to (re)initialize the counter as +# desired, or `gensym/counter` may be examined or altered directly +# to otherwise modify the name sequence generated by this function. +# +# Arguments: +# base_name = (in) (optional) The base name for the name generation. +# This argument defaults to `temp_` ifo +# omitted. +# +# Bugs: +# No checks are made for whether the names in question are indeed +# unassigned. +# +gensym := +proc(opt_base_name::string) +global `gensym/counter`; # in out +local base_name, tn; + +if (nargs >= 1) + then base_name := opt_base_name; + else base_name := '`temp_`'; +end if; + +# force an initialization if it hasn't been done before +if (not assigned(`gensym/counter`)) + then `gensym/init`(); +end if; + +tn := cat(base_name, `gensym/counter`); +`gensym/counter` := `gensym/counter` + 1; + +tn; return '%'; +end proc; + +############################################################################### + +# +# This procedure (re)initializes the gensym() global naming counter. +# +# Arguments: +# initial_counter = (in) (optional) +# The naming counter value to (re)initialize with. +# This argument defaults to 1 if omitted. +# +`gensym/init` := +proc(opt_initial_counter::integer) +global `gensym/counter`; # in out +local initial_counter; + +if (nargs >= 1) + then initial_counter := opt_initial_counter; + else initial_counter := 1; +end if; + +`gensym/counter` := initial_counter; +NULL; +end proc; + +############################################################################### +############################################################################### +############################################################################### + +# +# This function optionally saves a local variable to a global variable +# for later examination. This is useful for debugging purposes, since +# Maple doesn't have a proper debugger :-( . +# +# Arguments (explicit): +# fn_name = (in) The calling function or subsystem's name. +# expr = (in) The temporary result to save. +# name_postfix... = (in) (varargs) +# One or more arguments which this function +# concatenates to form the name postfix for +# the global variable. +# +# Arguments (implicit, as global variables) +# `saveit/level` = (in) (optional) If this name is assigned (as determined +# by the assigned() function), this +# function does the saving. Otherwise, +# this function is a nop. +# +saveit := +proc(n::integer, fn_name::string, label::string, expr::anything) +global `saveit/level`; +local save_name; + +if ( assigned(`saveit/level`) + and type(`saveit/level`, integer) + and (`saveit/level` >= n) ) + then + save_name := cat(fn_name,`/`,label); + lprint(printf(` --> \"%s\"`, save_name)); + assign( eval(save_name,1) = expr ); +end if; + +NULL; +end proc; |