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

1 files changed, 265 insertions, 11 deletions

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 $
 > 
 #
 # <<<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 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