From 8b734e46072a9614f8ffe12d3a10dcca558644a6 Mon Sep 17 00:00:00 2001 From: jthorn Date: Thu, 18 Jul 2002 17:37:11 +0000 Subject: directory to hold stuff that's no longer used but that someone might want to look at again git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@630 f88db872-0e4f-0410-b76b-b9085cfa78c5 --- archive/README | 6 +++ archive/ellipsoid.c | 28 ++++++++++ archive/ellipsoid.maple | 35 +++++++++++++ archive/ellipsoid.mm | 35 +++++++++++++ archive/ellipsoid.out | 38 ++++++++++++++ archive/setup_ellipsoid_initial_guess.cc | 87 ++++++++++++++++++++++++++++++++ 6 files changed, 229 insertions(+) create mode 100644 archive/README create mode 100644 archive/ellipsoid.c create mode 100644 archive/ellipsoid.maple create mode 100644 archive/ellipsoid.mm create mode 100644 archive/ellipsoid.out create mode 100644 archive/setup_ellipsoid_initial_guess.cc (limited to 'archive') diff --git a/archive/README b/archive/README new file mode 100644 index 0000000..6601932 --- /dev/null +++ b/archive/README @@ -0,0 +1,6 @@ +This directory contains stuff that's no longer used, but that I +thought might be useful for future reference, so I didn't want to +throw it away. + +Compared to the CVS attic, having a separate directory like this +seems more convenient -- you can 'cat, awk, and grep' more easily. diff --git a/archive/ellipsoid.c b/archive/ellipsoid.c new file mode 100644 index 0000000..49a6457 --- /dev/null +++ b/archive/ellipsoid.c @@ -0,0 +1,28 @@ +fp t1, t2, t3, t5, t6, t7, t9, t10, t12, t28; +fp t30, t33, t35, t36, t40, t42, t43, t48, t49, t52; +fp t55; + t1 = a*a; + t2 = b*b; + t3 = t1*t2; + t5 = t3*zcos*CW; + t6 = c*c; + t7 = t1*t6; + t9 = t7*ycos*BV; + t10 = t2*t6; + t12 = t10*xcos*AU; + t28 = xcos*xcos; + t30 = CW*CW; + t33 = BV*BV; + t35 = t10*t28; + t36 = ycos*ycos; + t40 = AU*AU; + t42 = t7*t36; + t43 = zcos*zcos; + t48 = t3*t43; + t49 = -2.0*t1*zcos*CW*ycos*BV-2.0*t2*zcos*CW*xcos*AU-2.0*t6*ycos*BV*xcos* +AU+t2*t28*t30+t6*t28*t33-t35+t1*t36*t30+t6*t36*t40-t42+t1*t43*t33+t2*t43*t40- +t48; + t52 = sqrt(-t3*t6*t49); + t55 = 1/(t35+t42+t48); + r_plus = (t5+t9+t12+t52)*t55; + r_minus = (t5+t9+t12-t52)*t55; diff --git a/archive/ellipsoid.maple b/archive/ellipsoid.maple new file mode 100644 index 0000000..0c485fa --- /dev/null +++ b/archive/ellipsoid.maple @@ -0,0 +1,35 @@ +# ellipsoid.maple -- compute equations for offset ellipsoid setup +# $Header$ + +# +# ellipsoid has center (A,B,C), radius (a,b,c) +# angular coordinate system has center (U,V,W) +# +# direction cosines wrt angular coordinate center are (alpha,beta,gamma) +# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos) +# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r) +# +# then the equation of the ellipsoid is +# (U+xcos*r - A)^2 (V+ycos*r - B)^2 (W+zcos*r - C)^2 +# ----------------- + ---------------- + ----------------- = 1 +# a^2 b^2 c^2 +# +# to solve this, we introduce intermediate variables +# AU = A - U +# BV = B - V +# CW = C - W +# +eqn := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1; + +read "../maple/util.mm"; +read "../maple/codegen2.mm"; + +[solve(eqn, r)]; +map(simplify, %); +[r_plus = %[1], r_minus = %[2]]; +solnlist := [codegen[optimize](%)]; + +ftruncate("ellipsoid.c"); +print_name_list_dcl(temps_in_eqnlist(solnlist, [r_plus,r_minus]), + "fp", "ellipsoid.c"); +codegen[C](solnlist, filename="ellipsoid.c"); diff --git a/archive/ellipsoid.mm b/archive/ellipsoid.mm new file mode 100644 index 0000000..d58b424 --- /dev/null +++ b/archive/ellipsoid.mm @@ -0,0 +1,35 @@ +# ellipsoid.maple -- compute equations for offset ellipsoid setup +# $Header: /usr/local/svn/cvs-repositories/numrelcvs/AEIThorns/AHFinderDirect/archive/ellipsoid.mm,v 1.1 2002-07-18 17:37:11 jthorn Exp $ + +# +# ellipsoid has center (A,B,C), radius (a,b,c) +# angular coordinate system has center (U,V,W) +# +# direction cosines wrt angular coordinate center are (alpha,beta,gamma) +# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos) +# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r) +# +# then the equation of the ellipsoid is +# (U+xcos*r - A)^2 (V+ycos*r - B)^2 (W+zcos*r - C)^2 +# ----------------- + ---------------- + ----------------- = 1 +# a^2 b^2 c^2 +# +# to solve this, we introduce intermediate variables +# AU = A - U +# BV = B - V +# CW = C - W +# +eqn := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1; + +read "../maple/util.mm"; +read "../maple/codegen2.mm"; + +[solve(eqn, r)]; +map(simplify, %); +[r_plus = %[1], r_minus = %[2]]; +solnlist := [codegen[optimize](%)]; + +ftruncate("ellipsoid.c"); +print_name_list_dcl(temps_in_eqnlist(solnlist, [r_plus,r_minus]), + "fp", "ellipsoid.c"); +codegen[C](solnlist, filename="ellipsoid.c"); diff --git a/archive/ellipsoid.out b/archive/ellipsoid.out new file mode 100644 index 0000000..4773c57 --- /dev/null +++ b/archive/ellipsoid.out @@ -0,0 +1,38 @@ + |\^/| Maple 7 (IBM INTEL LINUX) +._|\| |/|_. Copyright (c) 2001 by Waterloo Maple Inc. + \ MAPLE / All rights reserved. Maple is a registered trademark of + <____ ____> Waterloo Maple Inc. + | Type ? for help. +# ellipsoid.maple -- compute equations for offset ellipsoid setup +# $Header: /usr/local/svn/cvs-repositories/numrelcvs/AEIThorns/AHFinderDirect/archive/ellipsoid.out,v 1.1 2002-07-18 17:37:11 jthorn Exp $ +> +# +# ellipsoid has center (A,B,C), radius (a,b,c) +# angular coordinate system has center (U,V,W) +# +# direction cosines wrt angular coordinate center are (alpha,beta,gamma) +# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos) +# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r) +# +# then the equation of the ellipsoid is +# (U+xcos*r - A)^2 (V+ycos*r - B)^2 (W+zcos*r - C)^2 +# ----------------- + ---------------- + ----------------- = 1 +# a^2 b^2 c^2 +# +# to solve this, we introduce intermediate variables +# AU = A - U +# BV = B - V +# CW = C - W +# +> eqn := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1; + 2 2 2 + (xcos r - AU) (ycos r - BV) (zcos r - CW) + eqn := -------------- + -------------- + -------------- = 1 + 2 2 2 + a b c + +> +> read "../maple/util.mm"; +Error, unable to read `../maple/util.mm` +> quit +bytes used=129844, alloc=196572, time=0.03 diff --git a/archive/setup_ellipsoid_initial_guess.cc b/archive/setup_ellipsoid_initial_guess.cc new file mode 100644 index 0000000..f8ba391 --- /dev/null +++ b/archive/setup_ellipsoid_initial_guess.cc @@ -0,0 +1,87 @@ +// +// This function sets up an ellipsoid in the gridfn h, using the +// formulas in "ellipsoid.maple" and the Maple-generated C code in +// "ellipsoid.c": +// +// ellipsoid has center (A,B,C), radius (a,b,c) +// angular coordinate system has center (U,V,W) +// +// direction cosines wrt angular coordinate center are (xcos,ycos,zcos) +// i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r) +// +// then the equation of the ellipsoid is +// (U+xcos*r - A)^2 (V+ycos*r - B)^2 (W+zcos*r - C)^2 +// ----------------- + ---------------- + ----------------- = 1 +// a^2 b^2 c^2 +// +// to solve this, we introduce intermediate variables +// AU = A - U +// BV = B - V +// CW = C - W +// +namespace { +void setup_ellipsoid_initial_guess + (patch_system& ps, + fp global_center_x, fp global_center_y, fp global_center_z, + fp radius_x, fp radius_y, fp radius_z) +{ +CCTK_VInfo(CCTK_THORNSTRING, + " h = ellipsoid: global_center=(%g,%g,%g)", + global_center_x, global_center_y, global_center_z); +CCTK_VInfo(CCTK_THORNSTRING, + " radius=(%g,%g,%g)", + radius_x, radius_y, radius_z); + + for (int pn = 0 ; pn < ps.N_patches() ; ++pn) + { + patch& p = ps.ith_patch(pn); + + for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) + { + for (int isigma = p.min_isigma() ; + isigma <= p.max_isigma() ; + ++isigma) + { + const fp rho = p.rho_of_irho(irho); + const fp sigma = p.sigma_of_isigma(isigma); + fp xcos, ycos, zcos; + p.xyzcos_of_rho_sigma(rho,sigma, xcos,ycos,zcos); + + // set up variables used by Maple-generated code + const fp AU = global_center_x - ps.origin_x(); + const fp BV = global_center_y - ps.origin_y(); + const fp CW = global_center_z - ps.origin_z(); + const fp a = radius_x; + const fp b = radius_y; + const fp c = radius_z; + + // compute the solutions r_plus and r_minus + fp r_plus, r_minus; + #include "ellipsoid.c" + + // exactly one of the solutions (call it r) should be positive + fp r; + if ((r_plus > 0.0) && (r_minus < 0.0)) + then r = r_plus; + else if ((r_plus < 0.0) && (r_minus > 0.0)) + then r = r_minus; + else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING, + "\n" +" expected exactly one r>0 solution, got 0 or 2!\n" +" %s patch (irho,isigma)=(%d,%d) ==> (rho,sigma)=(%g,%g)\n" +" direction cosines (xcos,ycos,zcos)=(%g,%g,%g)\n" +" ==> r_plus=%g r_minus=%g\n" + , + p.name(), irho, isigma, + double(rho), double(sigma), + double(xcos), double(ycos), double(zcos), + double(r_plus), double(r_minus)); + /*NOTREACHED*/ + + // r = horizon radius at this grid point + p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) = r; + } + } + } +} + } -- cgit v1.2.3