From 197a844b61bb5f252f6164a00caf7ad712c1211e Mon Sep 17 00:00:00 2001 From: jthorn Date: Tue, 20 Aug 2002 16:53:41 +0000 Subject: add new files which will be (are) needed by Hermite interpolator git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@94 df1f8a13-aa1d-4dd4-9681-27ded5b42416 --- src/GeneralizedPolynomial-Uniform/common/1d.log | 276 +++++++++++++++- .../common/2d.cube.size6/coeff-I.dcl.c | 36 +++ .../common/2d.cube.size6/coeff-I.store.c | 36 +++ .../common/2d.cube.size6/coeff-dx.dcl.c | 36 +++ .../common/2d.cube.size6/coeff-dx.store.c | 36 +++ .../common/2d.cube.size6/coeff-dxx.dcl.c | 36 +++ .../common/2d.cube.size6/coeff-dxx.store.c | 36 +++ .../common/2d.cube.size6/coeff-dxy.dcl.c | 36 +++ .../common/2d.cube.size6/coeff-dxy.store.c | 36 +++ .../common/2d.cube.size6/coeff-dy.dcl.c | 36 +++ .../common/2d.cube.size6/coeff-dy.store.c | 36 +++ .../common/2d.cube.size6/coeff-dyy.dcl.c | 36 +++ .../common/2d.cube.size6/coeff-dyy.store.c | 36 +++ .../common/2d.cube.size6/data-var.assign.c | 36 +++ .../common/2d.cube.size6/data-var.dcl.c | 36 +++ .../common/2d.cube.size6/interp-I.compute.c | 37 +++ .../common/2d.cube.size6/interp-dx.compute.c | 37 +++ .../common/2d.cube.size6/interp-dxx.compute.c | 37 +++ .../common/2d.cube.size6/interp-dxy.compute.c | 37 +++ .../common/2d.cube.size6/interp-dy.compute.c | 37 +++ .../common/2d.cube.size6/interp-dyy.compute.c | 37 +++ src/GeneralizedPolynomial-Uniform/common/2d.log | 354 ++++++++++++++++++++- src/GeneralizedPolynomial-Uniform/common/2d.maple | 57 ++++ src/GeneralizedPolynomial-Uniform/common/3d.log | 298 +++++++++++++++-- .../common/cube_posns.maple | 4 + src/GeneralizedPolynomial-Uniform/common/makefile | 2 + 26 files changed, 1671 insertions(+), 46 deletions(-) create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.dcl.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.store.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.dcl.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.store.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.dcl.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.store.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.dcl.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.store.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.dcl.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.store.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.dcl.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.store.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.assign.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.dcl.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-I.compute.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dx.compute.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxx.compute.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxy.compute.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dy.compute.c create mode 100644 src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dyy.compute.c (limited to 'src') diff --git a/src/GeneralizedPolynomial-Uniform/common/1d.log b/src/GeneralizedPolynomial-Uniform/common/1d.log index c3e9168..bf3e9bf 100644 --- a/src/GeneralizedPolynomial-Uniform/common/1d.log +++ b/src/GeneralizedPolynomial-Uniform/common/1d.log @@ -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 +# $Id: interpolate.maple,v 1.7 2002/08/19 14:30:43 jthorn Exp $ > # # <<>> -# 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)}; @@ -535,6 +537,171 @@ local posn, data_eqns, coeff_eqns; return subs(coeff_eqns, eval(fn))(op(coord_list)) 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 + > ################################################################################ > @@ -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 1d interpolation schemes -# $Id: 1d.maple,v 1.4 2002/05/14 15:54:01 jthorn Exp $ +# $Id: 1d.maple,v 1.1.1.1 2002/08/18 15:05:38 jthorn Exp $ > ################################################################################ > @@ -1306,4 +1560,4 @@ data_var_list_1d_size7 := [ > ################################################################################ > quit -bytes used=716820, alloc=655240, time=0.16 +bytes used=797948, alloc=720764, time=0.12 diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.dcl.c new file mode 100644 index 0000000..0f93f17 --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.dcl.c @@ -0,0 +1,36 @@ +fp coeff_I_m2_m2, + coeff_I_m1_m2, + coeff_I_0_m2, + coeff_I_p1_m2, + coeff_I_p2_m2, + coeff_I_p3_m2, + coeff_I_m2_m1, + coeff_I_m1_m1, + coeff_I_0_m1, + coeff_I_p1_m1, + coeff_I_p2_m1, + coeff_I_p3_m1, + coeff_I_m2_0, + coeff_I_m1_0, + coeff_I_0_0, + coeff_I_p1_0, + coeff_I_p2_0, + coeff_I_p3_0, + coeff_I_m2_p1, + coeff_I_m1_p1, + coeff_I_0_p1, + coeff_I_p1_p1, + coeff_I_p2_p1, + coeff_I_p3_p1, + coeff_I_m2_p2, + coeff_I_m1_p2, + coeff_I_0_p2, + coeff_I_p1_p2, + coeff_I_p2_p2, + coeff_I_p3_p2, + coeff_I_m2_p3, + coeff_I_m1_p3, + coeff_I_0_p3, + coeff_I_p1_p3, + coeff_I_p2_p3, + coeff_I_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.store.c new file mode 100644 index 0000000..dfdbfda --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.store.c @@ -0,0 +1,36 @@ +COEFF(-2,-2) = coeff_I_m2_m2; +COEFF(-1,-2) = coeff_I_m1_m2; +COEFF(0,-2) = coeff_I_0_m2; +COEFF(1,-2) = coeff_I_p1_m2; +COEFF(2,-2) = coeff_I_p2_m2; +COEFF(3,-2) = coeff_I_p3_m2; +COEFF(-2,-1) = coeff_I_m2_m1; +COEFF(-1,-1) = coeff_I_m1_m1; +COEFF(0,-1) = coeff_I_0_m1; +COEFF(1,-1) = coeff_I_p1_m1; +COEFF(2,-1) = coeff_I_p2_m1; +COEFF(3,-1) = coeff_I_p3_m1; +COEFF(-2,0) = coeff_I_m2_0; +COEFF(-1,0) = coeff_I_m1_0; +COEFF(0,0) = coeff_I_0_0; +COEFF(1,0) = coeff_I_p1_0; +COEFF(2,0) = coeff_I_p2_0; +COEFF(3,0) = coeff_I_p3_0; +COEFF(-2,1) = coeff_I_m2_p1; +COEFF(-1,1) = coeff_I_m1_p1; +COEFF(0,1) = coeff_I_0_p1; +COEFF(1,1) = coeff_I_p1_p1; +COEFF(2,1) = coeff_I_p2_p1; +COEFF(3,1) = coeff_I_p3_p1; +COEFF(-2,2) = coeff_I_m2_p2; +COEFF(-1,2) = coeff_I_m1_p2; +COEFF(0,2) = coeff_I_0_p2; +COEFF(1,2) = coeff_I_p1_p2; +COEFF(2,2) = coeff_I_p2_p2; +COEFF(3,2) = coeff_I_p3_p2; +COEFF(-2,3) = coeff_I_m2_p3; +COEFF(-1,3) = coeff_I_m1_p3; +COEFF(0,3) = coeff_I_0_p3; +COEFF(1,3) = coeff_I_p1_p3; +COEFF(2,3) = coeff_I_p2_p3; +COEFF(3,3) = coeff_I_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.dcl.c new file mode 100644 index 0000000..bf2d3ef --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.dcl.c @@ -0,0 +1,36 @@ +fp coeff_dx_m2_m2, + coeff_dx_m1_m2, + coeff_dx_0_m2, + coeff_dx_p1_m2, + coeff_dx_p2_m2, + coeff_dx_p3_m2, + coeff_dx_m2_m1, + coeff_dx_m1_m1, + coeff_dx_0_m1, + coeff_dx_p1_m1, + coeff_dx_p2_m1, + coeff_dx_p3_m1, + coeff_dx_m2_0, + coeff_dx_m1_0, + coeff_dx_0_0, + coeff_dx_p1_0, + coeff_dx_p2_0, + coeff_dx_p3_0, + coeff_dx_m2_p1, + coeff_dx_m1_p1, + coeff_dx_0_p1, + coeff_dx_p1_p1, + coeff_dx_p2_p1, + coeff_dx_p3_p1, + coeff_dx_m2_p2, + coeff_dx_m1_p2, + coeff_dx_0_p2, + coeff_dx_p1_p2, + coeff_dx_p2_p2, + coeff_dx_p3_p2, + coeff_dx_m2_p3, + coeff_dx_m1_p3, + coeff_dx_0_p3, + coeff_dx_p1_p3, + coeff_dx_p2_p3, + coeff_dx_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.store.c new file mode 100644 index 0000000..233b5f3 --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.store.c @@ -0,0 +1,36 @@ +COEFF(-2,-2) = factor * coeff_dx_m2_m2; +COEFF(-1,-2) = factor * coeff_dx_m1_m2; +COEFF(0,-2) = factor * coeff_dx_0_m2; +COEFF(1,-2) = factor * coeff_dx_p1_m2; +COEFF(2,-2) = factor * coeff_dx_p2_m2; +COEFF(3,-2) = factor * coeff_dx_p3_m2; +COEFF(-2,-1) = factor * coeff_dx_m2_m1; +COEFF(-1,-1) = factor * coeff_dx_m1_m1; +COEFF(0,-1) = factor * coeff_dx_0_m1; +COEFF(1,-1) = factor * coeff_dx_p1_m1; +COEFF(2,-1) = factor * coeff_dx_p2_m1; +COEFF(3,-1) = factor * coeff_dx_p3_m1; +COEFF(-2,0) = factor * coeff_dx_m2_0; +COEFF(-1,0) = factor * coeff_dx_m1_0; +COEFF(0,0) = factor * coeff_dx_0_0; +COEFF(1,0) = factor * coeff_dx_p1_0; +COEFF(2,0) = factor * coeff_dx_p2_0; +COEFF(3,0) = factor * coeff_dx_p3_0; +COEFF(-2,1) = factor * coeff_dx_m2_p1; +COEFF(-1,1) = factor * coeff_dx_m1_p1; +COEFF(0,1) = factor * coeff_dx_0_p1; +COEFF(1,1) = factor * coeff_dx_p1_p1; +COEFF(2,1) = factor * coeff_dx_p2_p1; +COEFF(3,1) = factor * coeff_dx_p3_p1; +COEFF(-2,2) = factor * coeff_dx_m2_p2; +COEFF(-1,2) = factor * coeff_dx_m1_p2; +COEFF(0,2) = factor * coeff_dx_0_p2; +COEFF(1,2) = factor * coeff_dx_p1_p2; +COEFF(2,2) = factor * coeff_dx_p2_p2; +COEFF(3,2) = factor * coeff_dx_p3_p2; +COEFF(-2,3) = factor * coeff_dx_m2_p3; +COEFF(-1,3) = factor * coeff_dx_m1_p3; +COEFF(0,3) = factor * coeff_dx_0_p3; +COEFF(1,3) = factor * coeff_dx_p1_p3; +COEFF(2,3) = factor * coeff_dx_p2_p3; +COEFF(3,3) = factor * coeff_dx_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.dcl.c new file mode 100644 index 0000000..be3aa1c --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.dcl.c @@ -0,0 +1,36 @@ +fp coeff_dxx_m2_m2, + coeff_dxx_m1_m2, + coeff_dxx_0_m2, + coeff_dxx_p1_m2, + coeff_dxx_p2_m2, + coeff_dxx_p3_m2, + coeff_dxx_m2_m1, + coeff_dxx_m1_m1, + coeff_dxx_0_m1, + coeff_dxx_p1_m1, + coeff_dxx_p2_m1, + coeff_dxx_p3_m1, + coeff_dxx_m2_0, + coeff_dxx_m1_0, + coeff_dxx_0_0, + coeff_dxx_p1_0, + coeff_dxx_p2_0, + coeff_dxx_p3_0, + coeff_dxx_m2_p1, + coeff_dxx_m1_p1, + coeff_dxx_0_p1, + coeff_dxx_p1_p1, + coeff_dxx_p2_p1, + coeff_dxx_p3_p1, + coeff_dxx_m2_p2, + coeff_dxx_m1_p2, + coeff_dxx_0_p2, + coeff_dxx_p1_p2, + coeff_dxx_p2_p2, + coeff_dxx_p3_p2, + coeff_dxx_m2_p3, + coeff_dxx_m1_p3, + coeff_dxx_0_p3, + coeff_dxx_p1_p3, + coeff_dxx_p2_p3, + coeff_dxx_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.store.c new file mode 100644 index 0000000..08c3da8 --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.store.c @@ -0,0 +1,36 @@ +COEFF(-2,-2) = factor * coeff_dxx_m2_m2; +COEFF(-1,-2) = factor * coeff_dxx_m1_m2; +COEFF(0,-2) = factor * coeff_dxx_0_m2; +COEFF(1,-2) = factor * coeff_dxx_p1_m2; +COEFF(2,-2) = factor * coeff_dxx_p2_m2; +COEFF(3,-2) = factor * coeff_dxx_p3_m2; +COEFF(-2,-1) = factor * coeff_dxx_m2_m1; +COEFF(-1,-1) = factor * coeff_dxx_m1_m1; +COEFF(0,-1) = factor * coeff_dxx_0_m1; +COEFF(1,-1) = factor * coeff_dxx_p1_m1; +COEFF(2,-1) = factor * coeff_dxx_p2_m1; +COEFF(3,-1) = factor * coeff_dxx_p3_m1; +COEFF(-2,0) = factor * coeff_dxx_m2_0; +COEFF(-1,0) = factor * coeff_dxx_m1_0; +COEFF(0,0) = factor * coeff_dxx_0_0; +COEFF(1,0) = factor * coeff_dxx_p1_0; +COEFF(2,0) = factor * coeff_dxx_p2_0; +COEFF(3,0) = factor * coeff_dxx_p3_0; +COEFF(-2,1) = factor * coeff_dxx_m2_p1; +COEFF(-1,1) = factor * coeff_dxx_m1_p1; +COEFF(0,1) = factor * coeff_dxx_0_p1; +COEFF(1,1) = factor * coeff_dxx_p1_p1; +COEFF(2,1) = factor * coeff_dxx_p2_p1; +COEFF(3,1) = factor * coeff_dxx_p3_p1; +COEFF(-2,2) = factor * coeff_dxx_m2_p2; +COEFF(-1,2) = factor * coeff_dxx_m1_p2; +COEFF(0,2) = factor * coeff_dxx_0_p2; +COEFF(1,2) = factor * coeff_dxx_p1_p2; +COEFF(2,2) = factor * coeff_dxx_p2_p2; +COEFF(3,2) = factor * coeff_dxx_p3_p2; +COEFF(-2,3) = factor * coeff_dxx_m2_p3; +COEFF(-1,3) = factor * coeff_dxx_m1_p3; +COEFF(0,3) = factor * coeff_dxx_0_p3; +COEFF(1,3) = factor * coeff_dxx_p1_p3; +COEFF(2,3) = factor * coeff_dxx_p2_p3; +COEFF(3,3) = factor * coeff_dxx_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.dcl.c new file mode 100644 index 0000000..6b7f784 --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.dcl.c @@ -0,0 +1,36 @@ +fp coeff_dxy_m2_m2, + coeff_dxy_m1_m2, + coeff_dxy_0_m2, + coeff_dxy_p1_m2, + coeff_dxy_p2_m2, + coeff_dxy_p3_m2, + coeff_dxy_m2_m1, + coeff_dxy_m1_m1, + coeff_dxy_0_m1, + coeff_dxy_p1_m1, + coeff_dxy_p2_m1, + coeff_dxy_p3_m1, + coeff_dxy_m2_0, + coeff_dxy_m1_0, + coeff_dxy_0_0, + coeff_dxy_p1_0, + coeff_dxy_p2_0, + coeff_dxy_p3_0, + coeff_dxy_m2_p1, + coeff_dxy_m1_p1, + coeff_dxy_0_p1, + coeff_dxy_p1_p1, + coeff_dxy_p2_p1, + coeff_dxy_p3_p1, + coeff_dxy_m2_p2, + coeff_dxy_m1_p2, + coeff_dxy_0_p2, + coeff_dxy_p1_p2, + coeff_dxy_p2_p2, + coeff_dxy_p3_p2, + coeff_dxy_m2_p3, + coeff_dxy_m1_p3, + coeff_dxy_0_p3, + coeff_dxy_p1_p3, + coeff_dxy_p2_p3, + coeff_dxy_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.store.c new file mode 100644 index 0000000..32313f6 --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.store.c @@ -0,0 +1,36 @@ +COEFF(-2,-2) = factor * coeff_dxy_m2_m2; +COEFF(-1,-2) = factor * coeff_dxy_m1_m2; +COEFF(0,-2) = factor * coeff_dxy_0_m2; +COEFF(1,-2) = factor * coeff_dxy_p1_m2; +COEFF(2,-2) = factor * coeff_dxy_p2_m2; +COEFF(3,-2) = factor * coeff_dxy_p3_m2; +COEFF(-2,-1) = factor * coeff_dxy_m2_m1; +COEFF(-1,-1) = factor * coeff_dxy_m1_m1; +COEFF(0,-1) = factor * coeff_dxy_0_m1; +COEFF(1,-1) = factor * coeff_dxy_p1_m1; +COEFF(2,-1) = factor * coeff_dxy_p2_m1; +COEFF(3,-1) = factor * coeff_dxy_p3_m1; +COEFF(-2,0) = factor * coeff_dxy_m2_0; +COEFF(-1,0) = factor * coeff_dxy_m1_0; +COEFF(0,0) = factor * coeff_dxy_0_0; +COEFF(1,0) = factor * coeff_dxy_p1_0; +COEFF(2,0) = factor * coeff_dxy_p2_0; +COEFF(3,0) = factor * coeff_dxy_p3_0; +COEFF(-2,1) = factor * coeff_dxy_m2_p1; +COEFF(-1,1) = factor * coeff_dxy_m1_p1; +COEFF(0,1) = factor * coeff_dxy_0_p1; +COEFF(1,1) = factor * coeff_dxy_p1_p1; +COEFF(2,1) = factor * coeff_dxy_p2_p1; +COEFF(3,1) = factor * coeff_dxy_p3_p1; +COEFF(-2,2) = factor * coeff_dxy_m2_p2; +COEFF(-1,2) = factor * coeff_dxy_m1_p2; +COEFF(0,2) = factor * coeff_dxy_0_p2; +COEFF(1,2) = factor * coeff_dxy_p1_p2; +COEFF(2,2) = factor * coeff_dxy_p2_p2; +COEFF(3,2) = factor * coeff_dxy_p3_p2; +COEFF(-2,3) = factor * coeff_dxy_m2_p3; +COEFF(-1,3) = factor * coeff_dxy_m1_p3; +COEFF(0,3) = factor * coeff_dxy_0_p3; +COEFF(1,3) = factor * coeff_dxy_p1_p3; +COEFF(2,3) = factor * coeff_dxy_p2_p3; +COEFF(3,3) = factor * coeff_dxy_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.dcl.c new file mode 100644 index 0000000..489d18d --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.dcl.c @@ -0,0 +1,36 @@ +fp coeff_dy_m2_m2, + coeff_dy_m1_m2, + coeff_dy_0_m2, + coeff_dy_p1_m2, + coeff_dy_p2_m2, + coeff_dy_p3_m2, + coeff_dy_m2_m1, + coeff_dy_m1_m1, + coeff_dy_0_m1, + coeff_dy_p1_m1, + coeff_dy_p2_m1, + coeff_dy_p3_m1, + coeff_dy_m2_0, + coeff_dy_m1_0, + coeff_dy_0_0, + coeff_dy_p1_0, + coeff_dy_p2_0, + coeff_dy_p3_0, + coeff_dy_m2_p1, + coeff_dy_m1_p1, + coeff_dy_0_p1, + coeff_dy_p1_p1, + coeff_dy_p2_p1, + coeff_dy_p3_p1, + coeff_dy_m2_p2, + coeff_dy_m1_p2, + coeff_dy_0_p2, + coeff_dy_p1_p2, + coeff_dy_p2_p2, + coeff_dy_p3_p2, + coeff_dy_m2_p3, + coeff_dy_m1_p3, + coeff_dy_0_p3, + coeff_dy_p1_p3, + coeff_dy_p2_p3, + coeff_dy_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.store.c new file mode 100644 index 0000000..4ab90c2 --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.store.c @@ -0,0 +1,36 @@ +COEFF(-2,-2) = factor * coeff_dy_m2_m2; +COEFF(-1,-2) = factor * coeff_dy_m1_m2; +COEFF(0,-2) = factor * coeff_dy_0_m2; +COEFF(1,-2) = factor * coeff_dy_p1_m2; +COEFF(2,-2) = factor * coeff_dy_p2_m2; +COEFF(3,-2) = factor * coeff_dy_p3_m2; +COEFF(-2,-1) = factor * coeff_dy_m2_m1; +COEFF(-1,-1) = factor * coeff_dy_m1_m1; +COEFF(0,-1) = factor * coeff_dy_0_m1; +COEFF(1,-1) = factor * coeff_dy_p1_m1; +COEFF(2,-1) = factor * coeff_dy_p2_m1; +COEFF(3,-1) = factor * coeff_dy_p3_m1; +COEFF(-2,0) = factor * coeff_dy_m2_0; +COEFF(-1,0) = factor * coeff_dy_m1_0; +COEFF(0,0) = factor * coeff_dy_0_0; +COEFF(1,0) = factor * coeff_dy_p1_0; +COEFF(2,0) = factor * coeff_dy_p2_0; +COEFF(3,0) = factor * coeff_dy_p3_0; +COEFF(-2,1) = factor * coeff_dy_m2_p1; +COEFF(-1,1) = factor * coeff_dy_m1_p1; +COEFF(0,1) = factor * coeff_dy_0_p1; +COEFF(1,1) = factor * coeff_dy_p1_p1; +COEFF(2,1) = factor * coeff_dy_p2_p1; +COEFF(3,1) = factor * coeff_dy_p3_p1; +COEFF(-2,2) = factor * coeff_dy_m2_p2; +COEFF(-1,2) = factor * coeff_dy_m1_p2; +COEFF(0,2) = factor * coeff_dy_0_p2; +COEFF(1,2) = factor * coeff_dy_p1_p2; +COEFF(2,2) = factor * coeff_dy_p2_p2; +COEFF(3,2) = factor * coeff_dy_p3_p2; +COEFF(-2,3) = factor * coeff_dy_m2_p3; +COEFF(-1,3) = factor * coeff_dy_m1_p3; +COEFF(0,3) = factor * coeff_dy_0_p3; +COEFF(1,3) = factor * coeff_dy_p1_p3; +COEFF(2,3) = factor * coeff_dy_p2_p3; +COEFF(3,3) = factor * coeff_dy_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.dcl.c new file mode 100644 index 0000000..056404a --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.dcl.c @@ -0,0 +1,36 @@ +fp coeff_dyy_m2_m2, + coeff_dyy_m1_m2, + coeff_dyy_0_m2, + coeff_dyy_p1_m2, + coeff_dyy_p2_m2, + coeff_dyy_p3_m2, + coeff_dyy_m2_m1, + coeff_dyy_m1_m1, + coeff_dyy_0_m1, + coeff_dyy_p1_m1, + coeff_dyy_p2_m1, + coeff_dyy_p3_m1, + coeff_dyy_m2_0, + coeff_dyy_m1_0, + coeff_dyy_0_0, + coeff_dyy_p1_0, + coeff_dyy_p2_0, + coeff_dyy_p3_0, + coeff_dyy_m2_p1, + coeff_dyy_m1_p1, + coeff_dyy_0_p1, + coeff_dyy_p1_p1, + coeff_dyy_p2_p1, + coeff_dyy_p3_p1, + coeff_dyy_m2_p2, + coeff_dyy_m1_p2, + coeff_dyy_0_p2, + coeff_dyy_p1_p2, + coeff_dyy_p2_p2, + coeff_dyy_p3_p2, + coeff_dyy_m2_p3, + coeff_dyy_m1_p3, + coeff_dyy_0_p3, + coeff_dyy_p1_p3, + coeff_dyy_p2_p3, + coeff_dyy_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.store.c new file mode 100644 index 0000000..1a77da9 --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.store.c @@ -0,0 +1,36 @@ +COEFF(-2,-2) = factor * coeff_dyy_m2_m2; +COEFF(-1,-2) = factor * coeff_dyy_m1_m2; +COEFF(0,-2) = factor * coeff_dyy_0_m2; +COEFF(1,-2) = factor * coeff_dyy_p1_m2; +COEFF(2,-2) = factor * coeff_dyy_p2_m2; +COEFF(3,-2) = factor * coeff_dyy_p3_m2; +COEFF(-2,-1) = factor * coeff_dyy_m2_m1; +COEFF(-1,-1) = factor * coeff_dyy_m1_m1; +COEFF(0,-1) = factor * coeff_dyy_0_m1; +COEFF(1,-1) = factor * coeff_dyy_p1_m1; +COEFF(2,-1) = factor * coeff_dyy_p2_m1; +COEFF(3,-1) = factor * coeff_dyy_p3_m1; +COEFF(-2,0) = factor * coeff_dyy_m2_0; +COEFF(-1,0) = factor * coeff_dyy_m1_0; +COEFF(0,0) = factor * coeff_dyy_0_0; +COEFF(1,0) = factor * coeff_dyy_p1_0; +COEFF(2,0) = factor * coeff_dyy_p2_0; +COEFF(3,0) = factor * coeff_dyy_p3_0; +COEFF(-2,1) = factor * coeff_dyy_m2_p1; +COEFF(-1,1) = factor * coeff_dyy_m1_p1; +COEFF(0,1) = factor * coeff_dyy_0_p1; +COEFF(1,1) = factor * coeff_dyy_p1_p1; +COEFF(2,1) = factor * coeff_dyy_p2_p1; +COEFF(3,1) = factor * coeff_dyy_p3_p1; +COEFF(-2,2) = factor * coeff_dyy_m2_p2; +COEFF(-1,2) = factor * coeff_dyy_m1_p2; +COEFF(0,2) = factor * coeff_dyy_0_p2; +COEFF(1,2) = factor * coeff_dyy_p1_p2; +COEFF(2,2) = factor * coeff_dyy_p2_p2; +COEFF(3,2) = factor * coeff_dyy_p3_p2; +COEFF(-2,3) = factor * coeff_dyy_m2_p3; +COEFF(-1,3) = factor * coeff_dyy_m1_p3; +COEFF(0,3) = factor * coeff_dyy_0_p3; +COEFF(1,3) = factor * coeff_dyy_p1_p3; +COEFF(2,3) = factor * coeff_dyy_p2_p3; +COEFF(3,3) = factor * coeff_dyy_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.assign.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.assign.c new file mode 100644 index 0000000..0172f4a --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.assign.c @@ -0,0 +1,36 @@ +data_m2_m2 = DATA(-2,-2); +data_m1_m2 = DATA(-1,-2); +data_0_m2 = DATA(0,-2); +data_p1_m2 = DATA(1,-2); +data_p2_m2 = DATA(2,-2); +data_p3_m2 = DATA(3,-2); +data_m2_m1 = DATA(-2,-1); +data_m1_m1 = DATA(-1,-1); +data_0_m1 = DATA(0,-1); +data_p1_m1 = DATA(1,-1); +data_p2_m1 = DATA(2,-1); +data_p3_m1 = DATA(3,-1); +data_m2_0 = DATA(-2,0); +data_m1_0 = DATA(-1,0); +data_0_0 = DATA(0,0); +data_p1_0 = DATA(1,0); +data_p2_0 = DATA(2,0); +data_p3_0 = DATA(3,0); +data_m2_p1 = DATA(-2,1); +data_m1_p1 = DATA(-1,1); +data_0_p1 = DATA(0,1); +data_p1_p1 = DATA(1,1); +data_p2_p1 = DATA(2,1); +data_p3_p1 = DATA(3,1); +data_m2_p2 = DATA(-2,2); +data_m1_p2 = DATA(-1,2); +data_0_p2 = DATA(0,2); +data_p1_p2 = DATA(1,2); +data_p2_p2 = DATA(2,2); +data_p3_p2 = DATA(3,2); +data_m2_p3 = DATA(-2,3); +data_m1_p3 = DATA(-1,3); +data_0_p3 = DATA(0,3); +data_p1_p3 = DATA(1,3); +data_p2_p3 = DATA(2,3); +data_p3_p3 = DATA(3,3); diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.dcl.c new file mode 100644 index 0000000..342b513 --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.dcl.c @@ -0,0 +1,36 @@ +fp 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; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-I.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-I.compute.c new file mode 100644 index 0000000..90eaa22 --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-I.compute.c @@ -0,0 +1,37 @@ +result = + coeff_I_m2_m2*data_m2_m2 + + coeff_I_m1_m2*data_m1_m2 + + coeff_I_0_m2*data_0_m2 + + coeff_I_p1_m2*data_p1_m2 + + coeff_I_p2_m2*data_p2_m2 + + coeff_I_p3_m2*data_p3_m2 + + coeff_I_m2_m1*data_m2_m1 + + coeff_I_m1_m1*data_m1_m1 + + coeff_I_0_m1*data_0_m1 + + coeff_I_p1_m1*data_p1_m1 + + coeff_I_p2_m1*data_p2_m1 + + coeff_I_p3_m1*data_p3_m1 + + coeff_I_m2_0*data_m2_0 + + coeff_I_m1_0*data_m1_0 + + coeff_I_0_0*data_0_0 + + coeff_I_p1_0*data_p1_0 + + coeff_I_p2_0*data_p2_0 + + coeff_I_p3_0*data_p3_0 + + coeff_I_m2_p1*data_m2_p1 + + coeff_I_m1_p1*data_m1_p1 + + coeff_I_0_p1*data_0_p1 + + coeff_I_p1_p1*data_p1_p1 + + coeff_I_p2_p1*data_p2_p1 + + coeff_I_p3_p1*data_p3_p1 + + coeff_I_m2_p2*data_m2_p2 + + coeff_I_m1_p2*data_m1_p2 + + coeff_I_0_p2*data_0_p2 + + coeff_I_p1_p2*data_p1_p2 + + coeff_I_p2_p2*data_p2_p2 + + coeff_I_p3_p2*data_p3_p2 + + coeff_I_m2_p3*data_m2_p3 + + coeff_I_m1_p3*data_m1_p3 + + coeff_I_0_p3*data_0_p3 + + coeff_I_p1_p3*data_p1_p3 + + coeff_I_p2_p3*data_p2_p3 + + coeff_I_p3_p3*data_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dx.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dx.compute.c new file mode 100644 index 0000000..b67994e --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dx.compute.c @@ -0,0 +1,37 @@ +result = + coeff_dx_m2_m2*data_m2_m2 + + coeff_dx_m1_m2*data_m1_m2 + + coeff_dx_0_m2*data_0_m2 + + coeff_dx_p1_m2*data_p1_m2 + + coeff_dx_p2_m2*data_p2_m2 + + coeff_dx_p3_m2*data_p3_m2 + + coeff_dx_m2_m1*data_m2_m1 + + coeff_dx_m1_m1*data_m1_m1 + + coeff_dx_0_m1*data_0_m1 + + coeff_dx_p1_m1*data_p1_m1 + + coeff_dx_p2_m1*data_p2_m1 + + coeff_dx_p3_m1*data_p3_m1 + + coeff_dx_m2_0*data_m2_0 + + coeff_dx_m1_0*data_m1_0 + + coeff_dx_0_0*data_0_0 + + coeff_dx_p1_0*data_p1_0 + + coeff_dx_p2_0*data_p2_0 + + coeff_dx_p3_0*data_p3_0 + + coeff_dx_m2_p1*data_m2_p1 + + coeff_dx_m1_p1*data_m1_p1 + + coeff_dx_0_p1*data_0_p1 + + coeff_dx_p1_p1*data_p1_p1 + + coeff_dx_p2_p1*data_p2_p1 + + coeff_dx_p3_p1*data_p3_p1 + + coeff_dx_m2_p2*data_m2_p2 + + coeff_dx_m1_p2*data_m1_p2 + + coeff_dx_0_p2*data_0_p2 + + coeff_dx_p1_p2*data_p1_p2 + + coeff_dx_p2_p2*data_p2_p2 + + coeff_dx_p3_p2*data_p3_p2 + + coeff_dx_m2_p3*data_m2_p3 + + coeff_dx_m1_p3*data_m1_p3 + + coeff_dx_0_p3*data_0_p3 + + coeff_dx_p1_p3*data_p1_p3 + + coeff_dx_p2_p3*data_p2_p3 + + coeff_dx_p3_p3*data_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxx.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxx.compute.c new file mode 100644 index 0000000..cef82ff --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxx.compute.c @@ -0,0 +1,37 @@ +result = + coeff_dxx_m2_m2*data_m2_m2 + + coeff_dxx_m1_m2*data_m1_m2 + + coeff_dxx_0_m2*data_0_m2 + + coeff_dxx_p1_m2*data_p1_m2 + + coeff_dxx_p2_m2*data_p2_m2 + + coeff_dxx_p3_m2*data_p3_m2 + + coeff_dxx_m2_m1*data_m2_m1 + + coeff_dxx_m1_m1*data_m1_m1 + + coeff_dxx_0_m1*data_0_m1 + + coeff_dxx_p1_m1*data_p1_m1 + + coeff_dxx_p2_m1*data_p2_m1 + + coeff_dxx_p3_m1*data_p3_m1 + + coeff_dxx_m2_0*data_m2_0 + + coeff_dxx_m1_0*data_m1_0 + + coeff_dxx_0_0*data_0_0 + + coeff_dxx_p1_0*data_p1_0 + + coeff_dxx_p2_0*data_p2_0 + + coeff_dxx_p3_0*data_p3_0 + + coeff_dxx_m2_p1*data_m2_p1 + + coeff_dxx_m1_p1*data_m1_p1 + + coeff_dxx_0_p1*data_0_p1 + + coeff_dxx_p1_p1*data_p1_p1 + + coeff_dxx_p2_p1*data_p2_p1 + + coeff_dxx_p3_p1*data_p3_p1 + + coeff_dxx_m2_p2*data_m2_p2 + + coeff_dxx_m1_p2*data_m1_p2 + + coeff_dxx_0_p2*data_0_p2 + + coeff_dxx_p1_p2*data_p1_p2 + + coeff_dxx_p2_p2*data_p2_p2 + + coeff_dxx_p3_p2*data_p3_p2 + + coeff_dxx_m2_p3*data_m2_p3 + + coeff_dxx_m1_p3*data_m1_p3 + + coeff_dxx_0_p3*data_0_p3 + + coeff_dxx_p1_p3*data_p1_p3 + + coeff_dxx_p2_p3*data_p2_p3 + + coeff_dxx_p3_p3*data_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxy.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxy.compute.c new file mode 100644 index 0000000..7c66027 --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxy.compute.c @@ -0,0 +1,37 @@ +result = + coeff_dxy_m2_m2*data_m2_m2 + + coeff_dxy_m1_m2*data_m1_m2 + + coeff_dxy_0_m2*data_0_m2 + + coeff_dxy_p1_m2*data_p1_m2 + + coeff_dxy_p2_m2*data_p2_m2 + + coeff_dxy_p3_m2*data_p3_m2 + + coeff_dxy_m2_m1*data_m2_m1 + + coeff_dxy_m1_m1*data_m1_m1 + + coeff_dxy_0_m1*data_0_m1 + + coeff_dxy_p1_m1*data_p1_m1 + + coeff_dxy_p2_m1*data_p2_m1 + + coeff_dxy_p3_m1*data_p3_m1 + + coeff_dxy_m2_0*data_m2_0 + + coeff_dxy_m1_0*data_m1_0 + + coeff_dxy_0_0*data_0_0 + + coeff_dxy_p1_0*data_p1_0 + + coeff_dxy_p2_0*data_p2_0 + + coeff_dxy_p3_0*data_p3_0 + + coeff_dxy_m2_p1*data_m2_p1 + + coeff_dxy_m1_p1*data_m1_p1 + + coeff_dxy_0_p1*data_0_p1 + + coeff_dxy_p1_p1*data_p1_p1 + + coeff_dxy_p2_p1*data_p2_p1 + + coeff_dxy_p3_p1*data_p3_p1 + + coeff_dxy_m2_p2*data_m2_p2 + + coeff_dxy_m1_p2*data_m1_p2 + + coeff_dxy_0_p2*data_0_p2 + + coeff_dxy_p1_p2*data_p1_p2 + + coeff_dxy_p2_p2*data_p2_p2 + + coeff_dxy_p3_p2*data_p3_p2 + + coeff_dxy_m2_p3*data_m2_p3 + + coeff_dxy_m1_p3*data_m1_p3 + + coeff_dxy_0_p3*data_0_p3 + + coeff_dxy_p1_p3*data_p1_p3 + + coeff_dxy_p2_p3*data_p2_p3 + + coeff_dxy_p3_p3*data_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dy.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dy.compute.c new file mode 100644 index 0000000..ae07fcb --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dy.compute.c @@ -0,0 +1,37 @@ +result = + coeff_dy_m2_m2*data_m2_m2 + + coeff_dy_m1_m2*data_m1_m2 + + coeff_dy_0_m2*data_0_m2 + + coeff_dy_p1_m2*data_p1_m2 + + coeff_dy_p2_m2*data_p2_m2 + + coeff_dy_p3_m2*data_p3_m2 + + coeff_dy_m2_m1*data_m2_m1 + + coeff_dy_m1_m1*data_m1_m1 + + coeff_dy_0_m1*data_0_m1 + + coeff_dy_p1_m1*data_p1_m1 + + coeff_dy_p2_m1*data_p2_m1 + + coeff_dy_p3_m1*data_p3_m1 + + coeff_dy_m2_0*data_m2_0 + + coeff_dy_m1_0*data_m1_0 + + coeff_dy_0_0*data_0_0 + + coeff_dy_p1_0*data_p1_0 + + coeff_dy_p2_0*data_p2_0 + + coeff_dy_p3_0*data_p3_0 + + coeff_dy_m2_p1*data_m2_p1 + + coeff_dy_m1_p1*data_m1_p1 + + coeff_dy_0_p1*data_0_p1 + + coeff_dy_p1_p1*data_p1_p1 + + coeff_dy_p2_p1*data_p2_p1 + + coeff_dy_p3_p1*data_p3_p1 + + coeff_dy_m2_p2*data_m2_p2 + + coeff_dy_m1_p2*data_m1_p2 + + coeff_dy_0_p2*data_0_p2 + + coeff_dy_p1_p2*data_p1_p2 + + coeff_dy_p2_p2*data_p2_p2 + + coeff_dy_p3_p2*data_p3_p2 + + coeff_dy_m2_p3*data_m2_p3 + + coeff_dy_m1_p3*data_m1_p3 + + coeff_dy_0_p3*data_0_p3 + + coeff_dy_p1_p3*data_p1_p3 + + coeff_dy_p2_p3*data_p2_p3 + + coeff_dy_p3_p3*data_p3_p3; diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dyy.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dyy.compute.c new file mode 100644 index 0000000..dea1cd3 --- /dev/null +++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dyy.compute.c @@ -0,0 +1,37 @@ +result = + coeff_dyy_m2_m2*data_m2_m2 + + coeff_dyy_m1_m2*data_m1_m2 + + coeff_dyy_0_m2*data_0_m2 + + coeff_dyy_p1_m2*data_p1_m2 + + coeff_dyy_p2_m2*data_p2_m2 + + coeff_dyy_p3_m2*data_p3_m2 + + coeff_dyy_m2_m1*data_m2_m1 + + coeff_dyy_m1_m1*data_m1_m1 + + coeff_dyy_0_m1*data_0_m1 + + coeff_dyy_p1_m1*data_p1_m1 + + coeff_dyy_p2_m1*data_p2_m1 + + coeff_dyy_p3_m1*data_p3_m1 + + coeff_dyy_m2_0*data_m2_0 + + coeff_dyy_m1_0*data_m1_0 + + coeff_dyy_0_0*data_0_0 + + coeff_dyy_p1_0*data_p1_0 + + coeff_dyy_p2_0*data_p2_0 + + coeff_dyy_p3_0*data_p3_0 + + coeff_dyy_m2_p1*data_m2_p1 + + coeff_dyy_m1_p1*data_m1_p1 + + coeff_dyy_0_p1*data_0_p1 + + coeff_dyy_p1_p1*data_p1_p1 + + coeff_dyy_p2_p1*data_p2_p1 + + coeff_dyy_p3_p1*data_p3_p1 + + coeff_dyy_m2_p2*data_m2_p2 + + coeff_dyy_m1_p2*data_m1_p2 + + coeff_dyy_0_p2*data_0_p2 + + coeff_dyy_p1_p2*data_p1_p2 + + coeff_dyy_p2_p2*data_p2_p2 + + coeff_dyy_p3_p2*data_p3_p2 + + coeff_dyy_m2_p3*data_m2_p3 + + coeff_dyy_m1_p3*data_m1_p3 + + coeff_dyy_0_p3*data_0_p3 + + coeff_dyy_p1_p3*data_p1_p3 + + coeff_dyy_p2_p3*data_p2_p3 + + coeff_dyy_p3_p3*data_p3_p3; 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 $ > # # <<>> -# 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)}; @@ -535,6 +537,171 @@ local posn, data_eqns, coeff_eqns; return subs(coeff_eqns, eval(fn))(op(coord_list)) 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 + > ################################################################################ > @@ -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 diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.maple b/src/GeneralizedPolynomial-Uniform/common/2d.maple index 4a7b939..6efcc67 100644 --- a/src/GeneralizedPolynomial-Uniform/common/2d.maple +++ b/src/GeneralizedPolynomial-Uniform/common/2d.maple @@ -213,3 +213,60 @@ print_interp_cmpt__lc_of_data(posn_list_2d_size5, "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_"); + +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"); +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"); + +################################################################################ diff --git a/src/GeneralizedPolynomial-Uniform/common/3d.log b/src/GeneralizedPolynomial-Uniform/common/3d.log index 2fed690..12f76ca 100644 --- a/src/GeneralizedPolynomial-Uniform/common/3d.log +++ b/src/GeneralizedPolynomial-Uniform/common/3d.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 $ > # # <<>> -# 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)}; @@ -535,6 +537,171 @@ local posn, data_eqns, coeff_eqns; return subs(coeff_eqns, eval(fn))(op(coord_list)) 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 + > ################################################################################ > @@ -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 3d interpolation schemes -# $Id: 3d.maple,v 1.3 2002/05/14 15:54:01 jthorn Exp $ +# $Id: 3d.maple,v 1.1.1.1 2002/08/18 15:05:38 jthorn Exp $ > ################################################################################ > @@ -1168,12 +1422,12 @@ data_var_list_3d_size3 := ["data_m1_m1_m1", "data_0_m1_m1", "data_p1_m1_m1", > "3d.cube.size3/coeff-dyz.store.c"); > print_interp_coeff_var_store(posn_list_3d_size3, "factor", "coeff_dzz_", > "3d.cube.size3/coeff-dzz.store.c"); +bytes used=1000020, alloc=917336, time=0.16 > > print_name_list_dcl(map(coeff_name, posn_list_3d_size3, "coeff_I_"), "fp", > "3d.cube.size3/coeff-I.dcl.c"); > print_name_list_dcl(map(coeff_name, posn_list_3d_size3, "coeff_dx_"), "fp", > "3d.cube.size3/coeff-dx.dcl.c"); -bytes used=1000120, alloc=917336, time=0.18 > print_name_list_dcl(map(coeff_name, posn_list_3d_size3, "coeff_dy_"), "fp", > "3d.cube.size3/coeff-dy.dcl.c"); > print_name_list_dcl(map(coeff_name, posn_list_3d_size3, "coeff_dz_"), "fp", @@ -1229,6 +1483,7 @@ bytes used=1000120, alloc=917336, time=0.18 # > > data_var_list_3d_size4 := map(data_var_name, posn_list_3d_size4, "data_"); +bytes used=2000252, alloc=1179432, time=0.25 data_var_list_3d_size4 := ["data_m1_m1_m1", "data_0_m1_m1", "data_p1_m1_m1", "data_p2_m1_m1", "data_m1_0_m1", "data_0_0_m1", "data_p1_0_m1", @@ -1266,7 +1521,6 @@ data_var_list_3d_size4 := ["data_m1_m1_m1", "data_0_m1_m1", "data_p1_m1_m1", > "3d.cube.size4/data-var.dcl.c"); > print_data_var_assign(posn_list_3d_size4, "data_", > "3d.cube.size4/data-var.assign.c"); -bytes used=2000380, alloc=1179432, time=0.28 > > print_interp_coeff_var_store(posn_list_3d_size4, "", "coeff_I_", > "3d.cube.size4/coeff-I.store.c"); @@ -1293,9 +1547,9 @@ bytes used=2000380, alloc=1179432, time=0.28 > "3d.cube.size4/coeff-I.dcl.c"); > print_name_list_dcl(map(coeff_name, posn_list_3d_size4, "coeff_dx_"), "fp", > "3d.cube.size4/coeff-dx.dcl.c"); +bytes used=3000424, alloc=1244956, time=0.36 > print_name_list_dcl(map(coeff_name, posn_list_3d_size4, "coeff_dy_"), "fp", > "3d.cube.size4/coeff-dy.dcl.c"); -bytes used=3000632, alloc=1244956, time=0.40 > print_name_list_dcl(map(coeff_name, posn_list_3d_size4, "coeff_dz_"), "fp", > "3d.cube.size4/coeff-dz.dcl.c"); > print_name_list_dcl(map(coeff_name, posn_list_3d_size4, "coeff_dxx_"), "fp", @@ -1320,7 +1574,7 @@ bytes used=3000632, alloc=1244956, time=0.40 > print_interp_cmpt__lc_of_data(posn_list_3d_size4, > "result", "coeff_dy_", "data_", > "3d.cube.size4/interp-dy.compute.c"); -bytes used=4000828, alloc=1310480, time=0.49 +bytes used=4000728, alloc=1310480, time=0.46 > print_interp_cmpt__lc_of_data(posn_list_3d_size4, > "result", "coeff_dz_", "data_", > "3d.cube.size4/interp-dz.compute.c"); @@ -1342,6 +1596,7 @@ bytes used=4000828, alloc=1310480, time=0.49 > print_interp_cmpt__lc_of_data(posn_list_3d_size4, > "result", "coeff_dzz_", "data_", > "3d.cube.size4/interp-dzz.compute.c"); +bytes used=5000896, alloc=1376004, time=0.60 > ######################################## > @@ -1350,7 +1605,6 @@ bytes used=4000828, alloc=1310480, time=0.49 # > > data_var_list_3d_size5 := map(data_var_name, posn_list_3d_size5, "data_"); -bytes used=5001008, alloc=1310480, time=0.63 data_var_list_3d_size5 := ["data_m2_m2_m2", "data_m1_m2_m2", "data_0_m2_m2", "data_p1_m2_m2", "data_p2_m2_m2", "data_m2_m1_m2", "data_m1_m1_m2", @@ -1431,9 +1685,9 @@ data_var_list_3d_size5 := ["data_m2_m2_m2", "data_m1_m2_m2", "data_0_m2_m2", > "3d.cube.size5/coeff-dz.store.c"); > print_interp_coeff_var_store(posn_list_3d_size5, "factor", "coeff_dxx_", > "3d.cube.size5/coeff-dxx.store.c"); +bytes used=6001052, alloc=1376004, time=0.74 > print_interp_coeff_var_store(posn_list_3d_size5, "factor", "coeff_dxy_", > "3d.cube.size5/coeff-dxy.store.c"); -bytes used=6001252, alloc=1310480, time=0.77 > print_interp_coeff_var_store(posn_list_3d_size5, "factor", "coeff_dxz_", > "3d.cube.size5/coeff-dxz.store.c"); > print_interp_coeff_var_store(posn_list_3d_size5, "factor", "coeff_dyy_", @@ -1447,9 +1701,9 @@ bytes used=6001252, alloc=1310480, time=0.77 > "3d.cube.size5/coeff-I.dcl.c"); > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dx_"), "fp", > "3d.cube.size5/coeff-dx.dcl.c"); +bytes used=7001272, alloc=1376004, time=0.86 > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dy_"), "fp", > "3d.cube.size5/coeff-dy.dcl.c"); -bytes used=7001440, alloc=1376004, time=0.89 > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dz_"), "fp", > "3d.cube.size5/coeff-dz.dcl.c"); > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dxx_"), "fp", @@ -1462,7 +1716,7 @@ bytes used=7001440, alloc=1376004, time=0.89 > "3d.cube.size5/coeff-dyy.dcl.c"); > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dyz_"), "fp", > "3d.cube.size5/coeff-dyz.dcl.c"); -bytes used=8001776, alloc=1441528, time=0.98 +bytes used=8001444, alloc=1441528, time=0.98 > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dzz_"), "fp", > "3d.cube.size5/coeff-dzz.dcl.c"); > @@ -1475,10 +1729,10 @@ bytes used=8001776, alloc=1441528, time=0.98 > print_interp_cmpt__lc_of_data(posn_list_3d_size5, > "result", "coeff_dy_", "data_", > "3d.cube.size5/interp-dy.compute.c"); +bytes used=9001600, alloc=1441528, time=1.08 > print_interp_cmpt__lc_of_data(posn_list_3d_size5, > "result", "coeff_dz_", "data_", > "3d.cube.size5/interp-dz.compute.c"); -bytes used=9001948, alloc=1441528, time=1.09 > print_interp_cmpt__lc_of_data(posn_list_3d_size5, > "result", "coeff_dxx_", "data_", > "3d.cube.size5/interp-dxx.compute.c"); @@ -1488,10 +1742,10 @@ bytes used=9001948, alloc=1441528, time=1.09 > print_interp_cmpt__lc_of_data(posn_list_3d_size5, > "result", "coeff_dxz_", "data_", > "3d.cube.size5/interp-dxz.compute.c"); +bytes used=10001804, alloc=1507052, time=1.20 > print_interp_cmpt__lc_of_data(posn_list_3d_size5, > "result", "coeff_dyy_", "data_", > "3d.cube.size5/interp-dyy.compute.c"); -bytes used=10002200, alloc=1441528, time=1.21 > print_interp_cmpt__lc_of_data(posn_list_3d_size5, > "result", "coeff_dyz_", "data_", > "3d.cube.size5/interp-dyz.compute.c"); @@ -1501,4 +1755,4 @@ bytes used=10002200, alloc=1441528, time=1.21 > ################################################################################ > quit -bytes used=10722536, alloc=1441528, time=1.29 +bytes used=10794072, alloc=1507052, time=1.28 diff --git a/src/GeneralizedPolynomial-Uniform/common/cube_posns.maple b/src/GeneralizedPolynomial-Uniform/common/cube_posns.maple index 0d23e1c..13c134c 100644 --- a/src/GeneralizedPolynomial-Uniform/common/cube_posns.maple +++ b/src/GeneralizedPolynomial-Uniform/common/cube_posns.maple @@ -26,6 +26,8 @@ posn_list_2d_size4 := map(ListTools[Reverse], hypercube_points([-1,-1], [+2,+2])); posn_list_2d_size5 := map(ListTools[Reverse], hypercube_points([-2,-2], [+2,+2])); +posn_list_2d_size6 := map(ListTools[Reverse], + hypercube_points([-2,-2], [+3,+3])); ################################################################################ @@ -40,3 +42,5 @@ posn_list_3d_size4 := map(ListTools[Reverse], hypercube_points([-1,-1,-1], [+2,+2,+2])); posn_list_3d_size5 := map(ListTools[Reverse], hypercube_points([-2,-2,-2], [+2,+2,+2])); +posn_list_3d_size6 := map(ListTools[Reverse], + hypercube_points([-2,-2,-2], [+3,+3,+3])); diff --git a/src/GeneralizedPolynomial-Uniform/common/makefile b/src/GeneralizedPolynomial-Uniform/common/makefile index c34120b..77f2ead 100644 --- a/src/GeneralizedPolynomial-Uniform/common/makefile +++ b/src/GeneralizedPolynomial-Uniform/common/makefile @@ -52,6 +52,7 @@ dirs : 1d.dirs 2d.dirs 3d.dirs mkdir 2d.cube.size3 mkdir 2d.cube.size4 mkdir 2d.cube.size5 + mkdir 2d.cube.size6 .PHONY : 3d.dirs 3d.dirs: @@ -60,3 +61,4 @@ dirs : 1d.dirs 2d.dirs 3d.dirs mkdir 3d.cube.size3 mkdir 3d.cube.size4 mkdir 3d.cube.size5 + mkdir 3d.cube.size6 -- cgit v1.2.3