(* interpolatingList[n] is a function that creates n + 1 lists with list i representing the interpolating parameters for the ith basis function of degree at most n over a reference element of [-1, 1] with constraints on the value and derivative. The x values interpolated in each polynomial are {-1, 0, 1}. The function interpolates values depending on the conditions outlined in the background section of the paper. *)
interpolatingList[n_] :=
Module[{list, i, j}, 
list = ConstantArray[0, {n + 1, 3}] ; (* Define a list to hold the interpolant specifiers *)
For[i = 1, i <= n + 1, i++,
If[i == 1,
list[[i, 1]]= {-1, 1};
list[[i, 2]]=ConstantArray[0, n];
list[[i, 3]] = {1, 0};,
If[i == n+1,
list[[i, 1]]= {-1, 0};
list[[i, 2]]=ConstantArray[0, n];
list[[i, 3]] = {1, 1};,
list[[i, 1]]= {-1, 0};
list[[i, 2]]=ConstantArray[0, n];
list[[i, 2]][[i]] = 1;
list[[i, 3]] = {1, 0};
]
]
];
list
]

(* gen1DBasis[n, var] will generate the list of desired basis functions in var (e.g. x or y) of polynomial degree up to n. In the special case of n == 1, we do not interpolate at x = 0. *)
gen1DBasis[n_, var_]:= If[n == 1,{((1-var)/2),((1+var)/2)}, Array[Simplify[InterpolatingPolynomial[interpolatingList[n][[#]],var]]&, {n + 1}]];

(* genTable1D[n, var] is a helper function for genTable2D which makes sure that in the special case of n == 1, we can be sure that the number of basis functions line up properly when taking the 'product'. *) 
genTable1D[n_, var_] = If[n == 1, {((1-var)/2),0,((1+var)/2)}, gen1DBasis[n, var]];

(*genTable2D[p, q, n] generates a table with placeholder zeroes as described in the Implementation via Mathematica section of the paper. p is the degree of the basis functions in x, q is the degree of the basis functions in y, and n is the degree of the desired serendipity polynomials.
The format of the output is a table with each row representing an x polynomial which, in order by row, first interpolates at values x = -1 then interpolates at value x = 0, then with first, second, third, ..., (p - 1)st derivatives going from the top down,then the value at x = 1. Each column represents the same for the y polynomials going from left to right. At any point in the table, the entry is the product of the x and y polynomials which interpolate according to the row and column respectively. *)
genTable2D[p_, q_, n_] :=
Module[{table1,table2, pnew,qnew,i,j},
table1=
Which[
(* If p and q are both 1, adjoin zero rows and columns accounting for the linear case as discussed above. Due to the nature of this case, there are placeholder zeros for any functions which interpolate the value at at x = 0 or y = 0 *)
(p == 1&& q ==1),Join[Join[Table[genTable1D[p, x][[i]]*genTable1D[q, y][[j]], {i, 1, p +2}, {j, 1, q + 2}], ConstantArray[0,{p+2, Max[n-q-1, 0]}], 2], ConstantArray[0,{Max[n-p - 1, 0], n+1}]],
(* If p and q are the matching size for the target size of the table, generate the table with no placeholders. *)
(p ==n && q == n), Table[gen1DBasis[p, x][[i]]*gen1DBasis[q, y][[j]], {i, 1, p +1}, {j, 1, q + 1} ] , 
(* If p may be too small and q is 1, adjoin zero rows and columns accounting for the linear case as discussed above. *)
(1 <p ≤ n&& q  == 1),Join[Join[Table[genTable1D[p, x][[i]]*genTable1D[q, y][[j]], {i, 1, p +1}, {j, 1, q + 2}], ConstantArray[0, {p+1, n-q- 1}], 2], ConstantArray[0, {n-p, n+ 1}]],
(* If p is 1 and q may be too small, adjoin zero rows and columns accounting for the linear case as discussed above. *)
(p == 1&& 1 < q  ≤ n),Join[Join[Table[genTable1D[p, x][[i]]*genTable1D[q, y][[j]], {i, 1, p +2}, {j, 1, q + 1}], ConstantArray[0, {p+2, n-q}], 2], ConstantArray[0, {n-p - 1, n+1}]],
(* If at least one of p and q is too small, adjoin either zero rows or columns respectively to the end of the table. *)
(1 <  p ≤ n && 1  <  q ≤ n), Join[Join[Table[gen1DBasis[p, x][[i]]*gen1DBasis[q, y][[j]], {i, 1, p +1}, {j, 1, q + 1}], ConstantArray[0, {p+1, n-q}], 2], ConstantArray[0,{n-p, n+1}]] 
];
 (* Take into account the special linear case, changing p. *)
pnew=If[p==1,p+1,p];
 (* Take into account the special linear case, changing q. *) 
qnew=If[q==1,q+1,q];
(* Swap columns so that the zero columns are in the correct place i.e. in columns which interpolate derivatives, leaving the x = 1 value interpolating functions for the last column*)
table2=Join[Drop[table1,None,{qnew+1,qnew+1}],Take[table1,All,{qnew+1,qnew+1}],2];
 (* Swap rows so that the zero rows are in the correct place, in the same manner as the previous step. *)
Join[Drop[table2,{pnew+1}],Take[table2,{pnew+1}]]
]

(* The following give the basis functions for the tensor product spaces Q1, Q2, Q3, Q4, Q5, and Q6. *)

Q1 = genTable2D[1, 1, 1];

Q2 = genTable2D[2, 2, 2];

Q3 = genTable2D[3, 3, 3];

Q4 = genTable2D[4, 4, 4];

Q5 = genTable2D[5, 5, 5];

Q6 = genTable2D[6, 6, 6];

(* The following are the formulas which construct the basis functions for the serendipity spaces S1, S2, S3, S4, S5, and S6. *)

S1 = genTable2D[1, 1, 1];

S2 =Simplify[(genTable2D[1, 2, 2] ) + 
genTable2D[2, 1, 2 ] - genTable2D[1, 1, 2]] ;

S3 =Simplify[genTable2D[3, 1, 3]+
genTable2D[1, 3, 3]-
genTable2D[1, 1, 3]];

S4 =Simplify[genTable2D[1, 4, 4] +
genTable2D[4,1, 4]+genTable2D[2, 2, 4]+ 
(-1*genTable2D[1, 2, 4]) + 
(-1*genTable2D[2, 1, 4])];

S5=Simplify[genTable2D[1, 5, 5]+ 
genTable2D[5, 1, 5] + 
genTable2D[3, 2, 5] + 
genTable2D[2, 3, 5]+
(-1*genTable2D[1, 3, 5]) + 
(-1*genTable2D[3, 1, 5]) + 
(-1*genTable2D[2, 2, 5])];

S6 = Simplify[genTable2D[6, 1, 6] + 
genTable2D[1, 6, 6] + 
genTable2D[4, 2, 6] + 
genTable2D[2, 4, 6] + 
genTable2D[3, 3, 6] + 
(-1*genTable2D[4,1,6]) + 
(-1*genTable2D[1, 4, 6]) + 
(-1*genTable2D[2, 3, 6]) + 
(-1*genTable2D[3, 2, 6])];

(* The following is helper functions for genTable3D which will create a list of zeros of the desired length *) 
zeroLister3D[param1_, param2_, n_] := Table[Table[ConstantArray[0,n+1], {i, 1, param1}],{j, 1, param2}]

(*genTable3D[p, q, r, n] generates a table with placeholder zeroes. p is the degree of the basis functions in x, q is the degree of the basis functions in y, r is the degree of the basis functions in z, and n is desired serendipity degree.
The format of the output is a table with each row representing an x polynomial which, in order by row, first interpolates at values x = -1 then interpolates at value x = 0, then with first, second, third, ..., (p - 1)st derivatives going from the top down,then the value at x = 1. Each column represents the same for the y polynomials going from top down. The rows of inner lists represent the same for z polynomials. At any point in the table, the value output is the product of the x, y, and z polynomials corresponding to the interpolation properties of the row, column, and inner list row respectively. *)
genTable3D[p_, q_, r_, n_] :=
Module[{table1,table2,table3,pnew,qnew,rnew},
table1=
Which[
(* If p is 1, q is 1, and r is 1, adjoin zero rows, columns, and inner lists accounting for the linear case as discussed above. Again, there are extra placeholder zeros as a result of the construction. *)
(p == 1&& q  ==1 && r ==1 && n ==1),Join[Join[Join[Table[genTable1D[p, x][[i]]*genTable1D[q, y][[j]]*genTable1D[r, z][[k]], {i, 1, p +2}, {j, 1, q + 2}, {k, 1,r +2}], zeroLister3D[(q+2), (p+2), Max[n-r-2, 0]], 3], zeroLister3D[Max[n-2,0], p+2, n], 2], zeroLister3D[n + 1, Max[n-2,0], n]],
(* If p is 1, q is 1, and r is 1, adjoin zero rows, columns, and inner lists accounting for the linear case as discussed above. Again, there are extra placeholder zeros as a result of the construction. *)
(p == 1&& q  ==1 && r ==1),Join[Join[Join[Table[genTable1D[p, x][[i]]*genTable1D[q, y][[j]]*genTable1D[r, z][[k]], {i, 1, p +2}, {j, 1, q + 2}, {k, 1,r +2}], zeroLister3D[(q+2), (p+2), n-r-2], 3], zeroLister3D[n-2, p+2, n], 2], zeroLister3D[n + 1, n-2, n]],
(* If p, q, and r are the "correct" size, generate the table with no placeholder zeros. *)
(p ==n && q == n && r ==n),
 Table[gen1DBasis[p, x][[i]]*gen1DBasis[q, y][[j]]*gen1DBasis[r, z][[k]], {i, 1, p +1}, {j, 1, q + 1}, {k, 1, r + 1} ] ,
(* If p may be too small, q is 1, and r is 1, adjoin zero rows, columns, and inner lists accounting for the linear case as discussed above. *)
(1 <p ≤ n&& q  == 1 && r==1),Join[Join[Join[Table[genTable1D[p, x][[i]]*genTable1D[q, y][[j]]*genTable1D[r, z][[k]], {i, 1, p +1}, {j, 1, q + 2}, {k, 1, r + 2}], zeroLister3D[(q+2), (p+1), n-r-2], 3], zeroLister3D[n-2, p+1, n], 2], zeroLister3D[n+ 1, n-p, n]],
(* If p is 1, q may be too small, and r is 1, adjoin zero rows, columns, and inner lists accounting for the linear case as discussed above. *)
(p == 1&& 1 < q  ≤ n &&r==1),Join[Join[Join[Table[genTable1D[p, x][[i]]*genTable1D[q, y][[j]]*genTable1D[r, z][[k]], {i, 1, p +2}, {j, 1, q + 1}, {k, 1, r+2}],  zeroLister3D[(q+1), (p+2), n-r-2], 3], zeroLister3D[n-q, p+2, n], 2], zeroLister3D[n + 1, n-2, n]],
(* If p is 1, q is 1, and r may be too small, adjoin zero rows, columns, and inner lists accounting for the linear case as discussed above. *)
(p == 1&& q==1 &&1<r≤ n),Join[Join[Join[Table[genTable1D[p, x][[i]]*genTable1D[q, y][[j]]*genTable1D[r, z][[k]], {i, 1, p +2}, {j, 1, q + 2}, {k, 1, r+1}],  zeroLister3D[(q+2), (p+2), n-r-1], 3], zeroLister3D[n-2, p+2, n], 2], zeroLister3D[n + 1, n-2, n]],
(* If p may be too small, q may be too small, and r is 1, adjoin zero rows, columns, and inner lists accounting for the linear case as discussed above. *)
(1 <  p ≤ n && 1  <  q ≤ n && r==1), Join[Join[Join[Table[gen1DBasis[p, x][[i]]*gen1DBasis[q, y][[j]]*genTable1D[r, z][[k]], {i, 1, p +1}, {j, 1, q + 1}, {k, 1, r + 2}], zeroLister3D[(q+1), (p+1), n-r-2], 3], zeroLister3D[n-q, p+1, n], 2],zeroLister3D[n+ 1, n-p, n]],
(* If p may be too small, q is 1, and r may be too small, adjoin zero rows, columns, and inner lists accounting for the linear case as discussed above. *)
(1 <  p ≤ n && q==1 && 1<r≤ n), Join[Join[Join[Table[gen1DBasis[p, x][[i]]*genTable1D[q, y][[j]]*gen1DBasis[r, z][[k]], {i, 1, p +1}, {j, 1, q + 2}, {k, 1, r + 1}], zeroLister3D[(q+2), (p+1), n-r-1], 3], zeroLister3D[n-2, p+1, n], 2],zeroLister3D[n+ 1, n-p, n]],
(* If p is 1, q may be too small, and r may be too small, adjoin zero rows, columns, and inner lists accounting for the linear case as discussed above. *)
(p==1&& 1  <  q ≤ n &&1<r≤ n), Join[Join[Join[Table[genTable1D[p, x][[i]]*gen1DBasis[q, y][[j]]*gen1DBasis[r, z][[k]], {i, 1, p +2}, {j, 1, q + 1}, {k, 1, r + 1}], zeroLister3D[(q+1), (p+2), n-r-1], 3], zeroLister3D[n-q, p+2, n], 2],zeroLister3D[n+ 1, n-2, n]],
(* If p may be too small, q may be too small, and r may be too small, adjoin zero rows, columns, and inner lists accounting for the linear case as discussed above. *)
(1 <  p ≤ n && 1  <  q ≤ n && 1<r≤ n), Join[Join[Join[Table[gen1DBasis[p, x][[i]]*gen1DBasis[q, y][[j]]*gen1DBasis[r, z][[k]], {i, 1, p +1}, {j, 1, q + 1}, {k, 1, r + 1}], zeroLister3D[(q+1), (p+1), n-r-1], 3], zeroLister3D[n-q, p+1, n], 2],zeroLister3D[n+ 1, n-p, n]]
];
(* Take into account the special linear case, changing p. *)
pnew=If[p==1,p+1,p];
(* Take into account the special linear case, changing q. *)
qnew=If[q==1,q+1,q];
(* Take into account the special linear case, changing r. *)
rnew=If[r==1,r+1,r];
(* Swap inner lists so that the zero inner lists are in the correct place (see explanation for 2D case). *)
table2=Join[Drop[table1,None,None,{rnew+1,rnew+1}],Take[table1,All,All,{rnew+1,rnew+1}],3]; 
 (* Swap columns so that the zero columns are in the correct place. *)
table3=Join[Drop[table2,None,{qnew+1,qnew+1},None],Take[table2,All,{qnew+1,qnew+1},All],2];
(* Swap rows so that the zero rows are in the correct place. *)
Join[Drop[table3,{pnew+1,pnew+1},None,None],Take[table3,{pnew+1,pnew+1},All,All]] 
]

(* The following functions generates a table as a sum of tables given by permutations of a list as is used in the Floater-Gillette paper. *)

listPicker[alpha_, i_, n_] := genTable3D[Permutations[alpha][[i, 1]] , Permutations[alpha][[i, 2]], Permutations[alpha][[i, 3]], n]

q[alpha_, n_] :=Sum[listPicker[i, j, k], {i, {alpha}}, {k, {n}}, {j,Length[Permutations[alpha]]}]

(* The following give the basis functions for the tensor product spaces Q1, Q2, Q3, Q4, Q5, and Q6. *)

Q13D = genTable3D[1, 1, 1, 1];

Q23D = genTable3D[2, 2, 2, 2];

Q33D = genTable3D[3, 3, 3, 3];

Q43D = genTable3D[4, 4, 4, 4];

Q53D = genTable3D[5, 5, 5, 5];

Q63D = genTable3D[6, 6, 6, 6];

(*The following are formulas for the serendipity basis functions in three dimensions, S1, S2, S3, S4, and S5.*)

S13D = genTable3D[1,1,1,1];

S23D=Simplify[q[{1, 1, 2}, 2] -
2*q[{1, 1, 1}, 2]];

S33D=Simplify[q[{1, 1, 3},3]-
2*q[{1, 1, 1}, 3]];

S43D=Simplify[q[{1, 2, 2},4]+
(q[{1, 1, 4},4]+
(-2*q[{1, 1 ,2},4]))+
q[{1, 1, 1},4]];

S53D=Simplify[(q[{1, 2, 3} ,5]-
q[{1, 2, 2}, 5]) + 
(q[{1, 1, 5}, 5]+
(-2*q[{1, 1, 3}, 5])) + 
q[{1, 1, 1}, 5]];