--- /dev/null
+# Maple script to compute much of the data needed to implement the
+# family of Lagrange elements in 2d. Expects that the fields denoting
+# position and number of support points, etc are already set. Note that
+# we assume a bilinear mapping from the unit to the real cell.
+#
+# $Id$
+# Author: Wolfgang Bangerth, 1998
+
+
+
+ phi_polynom := array(0..n_functions-1);
+ grad_phi_polynom := array(0..n_functions-1);
+ local_mass_matrix := array(0..n_functions-1, 0..n_functions-1);
+
+ for i from 0 to n_functions-1 do
+ # note that the interp function wants vectors indexed from
+ # one and not from zero.
+ values := array(1..n_functions);
+ for j from 1 to n_functions do
+ values[j] := 0;
+ od;
+ values[i+1] := 1;
+
+ shifted_support_points := array (1..n_functions);
+ for j from 1 to n_functions do
+ shifted_support_points[j] := support_points[j-1];
+ od;
+
+ phi_polynom[i] := interp (shifted_support_points, values, xi);
+ grad_phi_polynom[i] := diff(phi_polynom[i], xi);
+ od;
+
+ phi:= proc(i,x) subs(xi=x, phi_polynom[i]); end;
+
+
+ points[0] := array(0..n_functions-1);
+ points[1] := array(0..n_functions-1);
+ for i from 0 to n_functions-1 do
+ points[0][i] := support_points[i]/2;
+ points[1][i] := support_points[i]/2+1/2;
+ od;
+
+ prolongation := array(0..1,0..n_functions-1, 0..n_functions-1);
+
+ for i from 0 to 1 do
+ for j from 0 to n_functions-1 do
+ for k from 0 to n_functions-1 do
+ prolongation[i,j,k] := phi(k, points[i][j]);
+ od;
+ od;
+ od;
+
+
+ # to get the restriction (interpolation) matrices, evaluate
+ # the basis functions on the child cells at the global
+ # interpolation points
+ child_phi[0] := proc(i, point)
+ if ((point<0) or (point>1/2)) then
+ 0:
+ else
+ phi(i,2*point):
+ fi:
+ end:
+ child_phi[1] := proc(i, point)
+ if ((point<1/2) or (point>1)) then
+ 0:
+ else
+ phi(i,2*point-1):
+ fi:
+ end:
+ restriction := array(0..1,0..n_functions-1, 0..n_functions-1);
+ for child from 0 to 1 do
+ for j from 0 to n_functions-1 do
+ for k from 0 to n_functions-1 do
+ restriction[child,j,k] := child_phi[child](k, support_points[j]):
+ od:
+ od:
+ od:
+
+
+ for i from 0 to n_functions-1 do
+ for j from 0 to n_functions-1 do
+ local_mass_matrix[i,j] := int(phi_polynom[i] * phi_polynom[j] * h,
+ xi=0..1);
+ od;
+ od;
+
--- /dev/null
+ n_functions := 4;
+
+ support_points := array(0..n_functions-1);
+ support_points[0] := 0;
+ support_points[1] := 1;
+ support_points[2] := 1/3;
+ support_points[3] := 2/3;
+
+
+ # do the real work
+ read "lagrange":
+
+
+ # write data to files
+ readlib(C);
+ C(phi_polynom, filename=cubic1d_shape_value);
+ C(grad_phi_polynom, filename=cubic1d_shape_grad);
+ C(prolongation, filename=cubic1d_prolongation);
+ C(restriction, filename=cubic1d_restriction);
+ C(local_mass_matrix, optimized, filename=cubic1d_massmatrix);
--- /dev/null
+ n_functions := 5;
+
+ support_points := array(0..n_functions-1);
+ support_points[0] := 0;
+ support_points[1] := 1;
+ support_points[2] := 1/4;
+ support_points[3] := 2/4;
+ support_points[4] := 3/4;
+
+
+ # do the real work
+ read "lagrange":
+
+
+ # write data to files
+ readlib(C);
+ C(phi_polynom, filename=cubic1d_shape_value);
+ C(grad_phi_polynom, filename=cubic1d_shape_grad);
+ C(prolongation, filename=cubic1d_prolongation);
+ C(restriction, filename=cubic1d_restriction);
+ C(local_mass_matrix, optimized, filename=cubic1d_massmatrix);
--- /dev/null
+# Use the following perl scripts to convert the output into a
+# suitable format:
+#
+# $Id$
+# Wolfgang Bangerth, 1998
+
+# concatenate lines belonging together
+perl -pi -e 's/([^;])\n/$1/g;' shape_value_1d
+perl -pi -e 's/([^;])\n/$1/g;' shape_grad_1d
+
+# give the programs a structure
+perl -pi -e 's/phi_polynom\[(\d)\] =/case $1: return/g;' *1d_shape_value
+perl -pi -e 's/grad_phi_polynom\[(\d)\] = (.*);/case $1: return Point<1>($2);/g;' *1d_shape_grad
+
+# use other indexing format for matrices
+perl -pi -e 's/\[(\d+)\]\[(\d)\]/($1,$2)/g;' *1d_massmatrix
+perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' *1d_prolongation
+perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' *1d_restriction
+
+# give temporaries a data type
+perl -pi -e 's/(t\d+) =/const double $1 =/g;' *1d_massmatrix
+
+# omit lines assigning zeroes to matrix elements, since zero is
+# already set and to save compilation time
+perl -pi -e 's/.*= 0.0;\n//g;' restriction_1d
+perl -pi -e 's/.*= 0.0;\n//g;' prolongation_1d
\ No newline at end of file