diff options
Diffstat (limited to 'src/GeneralizedPolynomial-Uniform/common/2d.log')
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/common/2d.log | 354 |
1 files changed, 341 insertions, 13 deletions
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.log b/src/GeneralizedPolynomial-Uniform/common/2d.log index 835c96e..909de57 100644 --- a/src/GeneralizedPolynomial-Uniform/common/2d.log +++ b/src/GeneralizedPolynomial-Uniform/common/2d.log @@ -4,7 +4,7 @@ <____ ____> Waterloo Maple Inc. | Type ? for help. # util.maple -- misc utility routines -# $Id: util.maple,v 1.3 2002/05/19 13:12:18 jthorn Exp $ +# $Header: /cactusdevcvs/CactusBase/LocalInterp/src/GeneralizedPolynomial-Uniform/util.maple,v 1.4 2002/08/20 16:46:06 jthorn Exp $ > # # fix_rationals - convert numbers to RATIONAL() calls @@ -444,12 +444,13 @@ ftruncate := proc(file_name::string) fopen(file_name, 'WRITE'); fclose(%); NULL end proc -# interpolate.maple -- compute generalized interpolation formulas/coefficients -# $Id: interpolate.maple,v 1.4 2002/05/14 15:52:50 jthorn Exp $ +# interpolate.maple -- compute interpolation formulas/coefficients +# $Header: /cactusdevcvs/CactusBase/LocalInterp/src/GeneralizedPolynomial-Uniform/interpolate.maple,v 1.8 2002/08/20 16:46:05 jthorn Exp $ > # # <<<representation of numbers, data values, etc>>> -# polynomial_interpolant - compute polynomial interpolant +# Lagrange_polynomial_interpolant - compute Lagrange polynomial interpolant +# Hermite_polynomial_interpolant - compute Hermite polynomial interpolant # coeff_as_lc_of_data - coefficients of ... (linear combination of data) # # print_coeff__lc_of_data - print C code to compute coefficients @@ -487,7 +488,8 @@ ftruncate := ################################################################################ > # -# This function computes a polynomial interpolant in any number of dimensions. +# This function computes a Lagrange polynomial interpolant in any +# number of dimensions. # # Arguments: # fn = The interpolation function. This should be a procedure in the @@ -507,7 +509,7 @@ ftruncate := # This function returns the interpolating polynomial, in the form of # an algebraic expression in the coordinates and the data values. # -> polynomial_interpolant := +> Lagrange_polynomial_interpolant := > proc( > fn::procedure, coeff_list::list(name), > coord_list::list(name), posn_list::list(list(numeric)) @@ -524,7 +526,7 @@ ftruncate := # interpolant as a polynomial in the coordinates > return subs(coeff_eqns, eval(fn))(op(coord_list)); > end proc; -polynomial_interpolant := proc(fn::procedure, coeff_list::list(name), +Lagrange_polynomial_interpolant := proc(fn::procedure, coeff_list::list(name), coord_list::list(name), posn_list::list(list(numeric))) local posn, data_eqns, coeff_eqns; data_eqns := {seq(fn(op(posn)) = 'DATA'(op(posn)), posn = posn_list)}; @@ -539,6 +541,171 @@ end proc ################################################################################ > # +# This function computes a Hermite polynomial interpolant in any +# number of dimensions. This is a polynomial which +# has values which match the given data DATA() at a specified set of +# points, and +# * has derivatives which match the specified finite-difference derivatives +# of the given data DATA() at a specified set of points +# +# For the derivative matching, we actually match all possible products +# of 1st derivatives, i.e. in 2-D we match dx, dy, and dxy, in 3-D we +# match dx, dy, dz, dxy, dxz, dyz, and dxyz, etc etc. +# +# Arguments: +# fn = The interpolation function. This should be a procedure in the +# coordinates, having the coefficients as global variables. For +# example, +# proc(x,y) +# + c03*y^3 + c13*x*y^3 + c23*x^2*y^3 + c33*x^3*y^3 +# + c02*y^2 + c12*x*y^2 + c22*x^2*y^2 + c32*x^3*y^2 +# + c01*y + c11*x*y + c21*x^2*y + c31*x^3*y +# + c00 + c10*x + c20*x^2 + c30*x^3 +# end proc; +# coeff_list = A set of the interpolation coefficients (coefficients in +# the interpolation function), for example +# [ +# c03, c13, c23, c33, +# c02, c12, c22, c32, +# c01, c11, c21, c31, +# c00, c10, c20, c30 +# ] +# coord_list = A list of the coordinates (independent variables in the +# interpolation function), for example [x,y]. +# deriv_coord_list = A list of lists of coordinates specifying which +# derivatives are computed by the deriv_proc_list[] +# procedures, for example +# [[x], [y], [x,y]] +# deriv_proc_list = A list of procedures for computing finite-difference +# derivatives. Each procedure should take N_dims integer +# arguments specifying an evaluation point, and return +# a suitable linear combination of the DATA() for the +# derivative at that point. For example +# example, +# [ +# proc(i::integer, j::integer) +# - 1/2*DATA(i-1,j) + 1/2*DATA(i+1,j) +# end proc +# , +# proc(i::integer, j::integer) +# - 1/2*DATA(i,j-1) + 1/2*DATA(i,j+1) +# end proc +# , +# proc(i::integer, j::integer) +# - 1/4*DATA(i-1,j+1) + 1/4*DATA(i+1,j+1) +# + 1/4*DATA(i-1,j-1) - 1/4*DATA(i+1,j-1) +# end proc +# ] +# fn_posn_list = A list of positions (each a list of numeric values) +# where the interpolant is to match the given data DATA(), +# for example +# [[0,0], [0,1], [1,0], [1,1]] +# deriv_posn_list = A list of positions (each a list of numeric values) +# where the interpolant is to match each possible product +# of 1st derivatives"1st" +# derivatives (actually all derivatives the given data DATA(), +# for example +# [[0,0], [0,1], [1,0], [1,1]] +# +# Results: +# This function returns the interpolating polynomial, in the form of +# an algebraic expression in the coordinates and the data values. +# +> Hermite_polynomial_interpolant := +> proc( +> fn::procedure, coeff_list::list(name), coord_list::list(name), +> deriv_coord_list::list(list(name)), fd_deriv_proc_list::list(procedure), +> fn_posn_list::list(list(numeric)), deriv_posn_list::list(list(numeric)) +> ) +> local fn_eset, +> fn_expr, deriv_eset; +> +# set of equations {fn(posn) = DATA(posn)} +> fn_eset := map( +> # return equation that fn(this point) = DATA(this point) +> proc(posn::list(integer)) +> fn(op(posn)) = 'DATA'(op(posn)); +> end proc +> , +> {op(fn_posn_list)} +> ); +> +# set of sets of equations +# { +# { deriv1(posn1) = fd_deriv1(posn1), +# deriv2(posn1) = fd_deriv2(posn1), +# ... }, +# { deriv1(posn2) = fd_deriv1(posn2), +# deriv2(posn2) = fd_deriv2(posn2), +# ... }, +# ... +# } +> fn_expr := fn(op(coord_list)); +> map( +> # return set of equations +> # {deriv(posn) = fd_deriv(posn)} +> # for this point +> proc(posn::list(integer)) +> { +> op( +> zip( +> # return equation that +> # deriv(posn) = fd_deriv(posn) +> # for this deriv and this point +> proc(deriv_coords::list(name), fd_deriv_proc::procedure) +> local fn_deriv_proc; +> fn_deriv_proc := unapply( diff(fn_expr,deriv_coords), +> op(coord_list) ); +> fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn)); +> end proc +> , +> [op(deriv_coord_list)] +> , +> [op(fd_deriv_proc_list)] +> ) +> ) +> } +> end proc +> , +> {op(deriv_posn_list)} +> ); +> +# set of equations {deriv-i(posn-j) = fd_deriv-i(posn-j)} +> deriv_eset := `union`(op(%)); +> +# solve equations for coeffs +> solve(fn_eset union deriv_eset, {op(coeff_list)}); +> +# interpolant as a polynomial in the coordinates +> return subs(%, eval(fn))(op(coord_list)); +> end proc; +Hermite_polynomial_interpolant := proc(fn::procedure, coeff_list::list(name), +coord_list::list(name), deriv_coord_list::list(list(name)), +fd_deriv_proc_list::list(procedure), fn_posn_list::list(list(numeric)), +deriv_posn_list::list(list(numeric))) +local fn_eset, fn_expr, deriv_eset; + fn_eset := map( + proc(posn::list(integer)) fn(op(posn)) = 'DATA'(op(posn)) end proc, + {op(fn_posn_list)}); + fn_expr := fn(op(coord_list)); + map(proc(posn::list(integer)) + {op(zip(proc(deriv_coords::list(name), fd_deriv_proc::procedure + ) + local fn_deriv_proc; + fn_deriv_proc := + unapply(diff(fn_expr, deriv_coords), op(coord_list)); + fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn)) + end proc, [op(deriv_coord_list)], [op(fd_deriv_proc_list)]))} + end proc, {op(deriv_posn_list)}); + deriv_eset := `union`(op(%)); + solve(fn_eset union deriv_eset, {op(coeff_list)}); + return subs(%, eval(fn))(op(coord_list)) +end proc + +> +################################################################################ +> +# # This function takes as input an interpolating polynomial, expresses # it as a linear combination of the data values, and returns the coefficeints # of that form. @@ -546,8 +713,7 @@ end proc # Arguments: # interpolant = The interpolating polynomial (an algebraic expression # in the coordinates and the data values). -# posn_list = The same list of positions as was used to compute the -# interpolating polynomial. +# posn_list = The same list of data positions used in the interpolant. # # Results: # This function returns the coefficients, as a list of equations of the @@ -915,7 +1081,7 @@ data_var_name := proc(posn::list(numeric), name_prefix::string) end proc # Maple code to compute lists of point positions in hypercube-shaped molecules -# $Id: $ +# $Id: cube_posns.maple,v 1.1.1.1 2002/08/18 15:05:38 jthorn Exp $ > ################################################################################ > @@ -974,6 +1140,18 @@ posn_list_2d_size5 := [[-2, -2], [-1, -2], [0, -2], [1, -2], [2, -2], [-2, -1], [1, 2], [2, 2]] +> posn_list_2d_size6 := map(ListTools[Reverse], +> hypercube_points([-2,-2], [+3,+3])); +posn_list_2d_size6 := [[-2, -2], [-1, -2], [0, -2], [1, -2], [2, -2], [3, -2], + + [-2, -1], [-1, -1], [0, -1], [1, -1], [2, -1], [3, -1], [-2, 0], [-1, 0], + + [0, 0], [1, 0], [2, 0], [3, 0], [-2, 1], [-1, 1], [0, 1], [1, 1], [2, 1], + + [3, 1], [-2, 2], [-1, 2], [0, 2], [1, 2], [2, 2], [3, 2], [-2, 3], [-1, 3], + + [0, 3], [1, 3], [2, 3], [3, 3]] + > ################################################################################ > @@ -1068,8 +1246,84 @@ posn_list_3d_size5 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2], [0, 2, 2], [1, 2, 2], [2, 2, 2]] +> posn_list_3d_size6 := map(ListTools[Reverse], +> hypercube_points([-2,-2,-2], [+3,+3,+3])); +posn_list_3d_size6 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2], + + [2, -2, -2], [3, -2, -2], [-2, -1, -2], [-1, -1, -2], [0, -1, -2], + + [1, -1, -2], [2, -1, -2], [3, -1, -2], [-2, 0, -2], [-1, 0, -2], [0, 0, -2], + + [1, 0, -2], [2, 0, -2], [3, 0, -2], [-2, 1, -2], [-1, 1, -2], [0, 1, -2], + + [1, 1, -2], [2, 1, -2], [3, 1, -2], [-2, 2, -2], [-1, 2, -2], [0, 2, -2], + + [1, 2, -2], [2, 2, -2], [3, 2, -2], [-2, 3, -2], [-1, 3, -2], [0, 3, -2], + + [1, 3, -2], [2, 3, -2], [3, 3, -2], [-2, -2, -1], [-1, -2, -1], [0, -2, -1], + + [1, -2, -1], [2, -2, -1], [3, -2, -1], [-2, -1, -1], [-1, -1, -1], + + [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [-2, 0, -1], + + [-1, 0, -1], [0, 0, -1], [1, 0, -1], [2, 0, -1], [3, 0, -1], [-2, 1, -1], + + [-1, 1, -1], [0, 1, -1], [1, 1, -1], [2, 1, -1], [3, 1, -1], [-2, 2, -1], + + [-1, 2, -1], [0, 2, -1], [1, 2, -1], [2, 2, -1], [3, 2, -1], [-2, 3, -1], + + [-1, 3, -1], [0, 3, -1], [1, 3, -1], [2, 3, -1], [3, 3, -1], [-2, -2, 0], + + [-1, -2, 0], [0, -2, 0], [1, -2, 0], [2, -2, 0], [3, -2, 0], [-2, -1, 0], + + [-1, -1, 0], [0, -1, 0], [1, -1, 0], [2, -1, 0], [3, -1, 0], [-2, 0, 0], + + [-1, 0, 0], [0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [-2, 1, 0], + + [-1, 1, 0], [0, 1, 0], [1, 1, 0], [2, 1, 0], [3, 1, 0], [-2, 2, 0], + + [-1, 2, 0], [0, 2, 0], [1, 2, 0], [2, 2, 0], [3, 2, 0], [-2, 3, 0], + + [-1, 3, 0], [0, 3, 0], [1, 3, 0], [2, 3, 0], [3, 3, 0], [-2, -2, 1], + + [-1, -2, 1], [0, -2, 1], [1, -2, 1], [2, -2, 1], [3, -2, 1], [-2, -1, 1], + + [-1, -1, 1], [0, -1, 1], [1, -1, 1], [2, -1, 1], [3, -1, 1], [-2, 0, 1], + + [-1, 0, 1], [0, 0, 1], [1, 0, 1], [2, 0, 1], [3, 0, 1], [-2, 1, 1], + + [-1, 1, 1], [0, 1, 1], [1, 1, 1], [2, 1, 1], [3, 1, 1], [-2, 2, 1], + + [-1, 2, 1], [0, 2, 1], [1, 2, 1], [2, 2, 1], [3, 2, 1], [-2, 3, 1], + + [-1, 3, 1], [0, 3, 1], [1, 3, 1], [2, 3, 1], [3, 3, 1], [-2, -2, 2], + + [-1, -2, 2], [0, -2, 2], [1, -2, 2], [2, -2, 2], [3, -2, 2], [-2, -1, 2], + + [-1, -1, 2], [0, -1, 2], [1, -1, 2], [2, -1, 2], [3, -1, 2], [-2, 0, 2], + + [-1, 0, 2], [0, 0, 2], [1, 0, 2], [2, 0, 2], [3, 0, 2], [-2, 1, 2], + + [-1, 1, 2], [0, 1, 2], [1, 1, 2], [2, 1, 2], [3, 1, 2], [-2, 2, 2], + + [-1, 2, 2], [0, 2, 2], [1, 2, 2], [2, 2, 2], [3, 2, 2], [-2, 3, 2], + + [-1, 3, 2], [0, 3, 2], [1, 3, 2], [2, 3, 2], [3, 3, 2], [-2, -2, 3], + + [-1, -2, 3], [0, -2, 3], [1, -2, 3], [2, -2, 3], [3, -2, 3], [-2, -1, 3], + + [-1, -1, 3], [0, -1, 3], [1, -1, 3], [2, -1, 3], [3, -1, 3], [-2, 0, 3], + + [-1, 0, 3], [0, 0, 3], [1, 0, 3], [2, 0, 3], [3, 0, 3], [-2, 1, 3], + + [-1, 1, 3], [0, 1, 3], [1, 1, 3], [2, 1, 3], [3, 1, 3], [-2, 2, 3], + + [-1, 2, 3], [0, 2, 3], [1, 2, 3], [2, 2, 3], [3, 2, 3], [-2, 3, 3], + + [-1, 3, 3], [0, 3, 3], [1, 3, 3], [2, 3, 3], [3, 3, 3]] + # Maple code to compute common coefficients for all 2d interpolation schemes -# $Id: 2d.maple,v 1.3 2002/05/14 15:54:01 jthorn Exp $ +# $Id: 2d.maple,v 1.1.1.1 2002/08/18 15:05:38 jthorn Exp $ > ################################################################################ > @@ -1218,6 +1472,7 @@ data_var_list_2d_size4 := ["data_m1_m1", "data_0_m1", "data_p1_m1", > "2d.cube.size4/coeff-dxy.dcl.c"); > print_name_list_dcl(map(coeff_name, posn_list_2d_size4, "coeff_dyy_"), "fp", > "2d.cube.size4/coeff-dyy.dcl.c"); +bytes used=1000212, alloc=917336, time=0.13 > > print_interp_cmpt__lc_of_data(posn_list_2d_size4, > "result", "coeff_I_", "data_", @@ -1228,7 +1483,6 @@ data_var_list_2d_size4 := ["data_m1_m1", "data_0_m1", "data_p1_m1", > print_interp_cmpt__lc_of_data(posn_list_2d_size4, > "result", "coeff_dy_", "data_", > "2d.cube.size4/interp-dy.compute.c"); -bytes used=1000024, alloc=917336, time=0.17 > print_interp_cmpt__lc_of_data(posn_list_2d_size4, > "result", "coeff_dxx_", "data_", > "2d.cube.size4/interp-dxx.compute.c"); @@ -1310,5 +1564,79 @@ data_var_list_2d_size5 := ["data_m2_m2", "data_m1_m2", "data_0_m2", > "2d.cube.size5/interp-dyy.compute.c"); > ################################################################################ +# +# generic stuff for 2d, cube, size=5 +# +> +> data_var_list_2d_size6 := map(data_var_name, posn_list_2d_size6, "data_"); +data_var_list_2d_size6 := ["data_m2_m2", "data_m1_m2", "data_0_m2", + + "data_p1_m2", "data_p2_m2", "data_p3_m2", "data_m2_m1", "data_m1_m1", + + "data_0_m1", "data_p1_m1", "data_p2_m1", "data_p3_m1", "data_m2_0", + + "data_m1_0", "data_0_0", "data_p1_0", "data_p2_0", "data_p3_0", + + "data_m2_p1", "data_m1_p1", "data_0_p1", "data_p1_p1", "data_p2_p1", + + "data_p3_p1", "data_m2_p2", "data_m1_p2", "data_0_p2", "data_p1_p2", + + "data_p2_p2", "data_p3_p2", "data_m2_p3", "data_m1_p3", "data_0_p3", + + "data_p1_p3", "data_p2_p3", "data_p3_p3"] + +> +> print_name_list_dcl(data_var_list_2d_size6, "fp", +> "2d.cube.size6/data-var.dcl.c"); +> print_data_var_assign(posn_list_2d_size6, "data_", +> "2d.cube.size6/data-var.assign.c"); +> +> print_interp_coeff_var_store(posn_list_2d_size6, "", "coeff_I_", +> "2d.cube.size6/coeff-I.store.c"); +> print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dx_", +> "2d.cube.size6/coeff-dx.store.c"); +> print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dy_", +> "2d.cube.size6/coeff-dy.store.c"); +> print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dxx_", +> "2d.cube.size6/coeff-dxx.store.c"); +bytes used=2000440, alloc=1179432, time=0.27 +> print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dxy_", +> "2d.cube.size6/coeff-dxy.store.c"); +> print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dyy_", +> "2d.cube.size6/coeff-dyy.store.c"); +> +> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_I_"), "fp", +> "2d.cube.size6/coeff-I.dcl.c"); +> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dx_"), "fp", +> "2d.cube.size6/coeff-dx.dcl.c"); +> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dy_"), "fp", +> "2d.cube.size6/coeff-dy.dcl.c"); +> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dxx_"), "fp", +> "2d.cube.size6/coeff-dxx.dcl.c"); +> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dxy_"), "fp", +> "2d.cube.size6/coeff-dxy.dcl.c"); +> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dyy_"), "fp", +> "2d.cube.size6/coeff-dyy.dcl.c"); +> +> print_interp_cmpt__lc_of_data(posn_list_2d_size6, +> "result", "coeff_I_", "data_", +> "2d.cube.size6/interp-I.compute.c"); +> print_interp_cmpt__lc_of_data(posn_list_2d_size6, +> "result", "coeff_dx_", "data_", +> "2d.cube.size6/interp-dx.compute.c"); +> print_interp_cmpt__lc_of_data(posn_list_2d_size6, +> "result", "coeff_dy_", "data_", +> "2d.cube.size6/interp-dy.compute.c"); +> print_interp_cmpt__lc_of_data(posn_list_2d_size6, +> "result", "coeff_dxx_", "data_", +> "2d.cube.size6/interp-dxx.compute.c"); +> print_interp_cmpt__lc_of_data(posn_list_2d_size6, +> "result", "coeff_dxy_", "data_", +> "2d.cube.size6/interp-dxy.compute.c"); +> print_interp_cmpt__lc_of_data(posn_list_2d_size6, +> "result", "coeff_dyy_", "data_", +> "2d.cube.size6/interp-dyy.compute.c"); +> +################################################################################ > quit -bytes used=1720848, alloc=917336, time=0.31 +bytes used=2688044, alloc=1179432, time=0.39 |