--- /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.
+#
+# $Id$
+# Author: Wolfgang Bangerth, 1998
+
+ phi_polynom := array(0..n_functions-1):
+ grad_phi_polynom := array(0..n_functions-1,0..1):
+ local_mass_matrix := array(0..n_functions-1, 0..n_functions-1):
+ prolongation := array(0..3,0..n_functions-1, 0..n_functions-1):
+ interface_constraints := array(0..2*(n_face_functions-2)+1-1,
+ 0..n_face_functions-1):
+ real_points := array(0..n_functions-1, 0..1);
+
+ print ("Computing basis functions"):
+ for i from 0 to n_functions-1 do
+ print (i):
+ values := array(1..n_functions):
+ for j from 1 to n_functions do
+ values[j] := 0:
+ od:
+ values[i+1] := 1:
+
+ equation_system := {}:
+ for j from 0 to n_functions-1 do
+ poly := subs(xi=support_points[j][1],
+ eta=support_points[j][2],
+ trial_function):
+ if (i=j) then
+ equation_system := equation_system union {poly = 1}:
+ else
+ equation_system := equation_system union {poly = 0}:
+ fi:
+ od:
+
+ phi_polynom[i] := subs(solve(equation_system), trial_function):
+ grad_phi_polynom[i,0] := diff(phi_polynom[i], xi):
+ grad_phi_polynom[i,1] := diff(phi_polynom[i], eta):
+ od:
+
+ phi:= proc(i,x,y) subs(xi=x, eta=y, phi_polynom[i]): end:
+
+
+ #points on children: let them be indexed one-based, as are
+ #the support_points
+ points[0] := array(0..n_functions-1, 1..2):
+ points[1] := array(0..n_functions-1, 1..2):
+ points[2] := array(0..n_functions-1, 1..2):
+ points[3] := array(0..n_functions-1, 1..2):
+ for i from 0 to n_functions-1 do
+ points[0][i,1] := support_points[i][1]/2:
+ points[0][i,2] := support_points[i][2]/2:
+
+ points[1][i,1] := support_points[i][1]/2+1/2:
+ points[1][i,2] := support_points[i][2]/2:
+
+ points[2][i,1] := support_points[i][1]/2+1/2:
+ points[2][i,2] := support_points[i][2]/2+1/2:
+
+ points[3][i,1] := support_points[i][1]/2:
+ points[3][i,2] := support_points[i][2]/2+1/2:
+ od:
+
+ print ("Computing prolongation matrices"):
+ for i from 0 to 3 do
+ print ("child", i):
+ 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,1], points[i][j,2]):
+ od:
+ od:
+ od:
+
+ print ("Computing restriction matrices"):
+ # to get the restriction (interpolation) matrices, evaluate
+ # the basis functions on the child cells at the global
+ # interpolation points
+ child_phi[0] := proc(i, x, y)
+ if ((x>1/2) or (y>1/2)) then
+ 0:
+ else
+ phi(i,2*x,2*y):
+ fi:
+ end:
+ child_phi[1] := proc(i, x, y)
+ if ((x<1/2) or (y>1/2)) then
+ 0:
+ else
+ phi(i,2*x-1,2*y):
+ fi:
+ end:
+ child_phi[2] := proc(i, x, y)
+ if ((x<1/2) or (y<1/2)) then
+ 0:
+ else
+ phi(i,2*x-1,2*y-1):
+ fi:
+ end:
+ child_phi[3] := proc(i, x, y)
+ if ((x>1/2) or (y<1/2)) then
+ 0:
+ else
+ phi(i,2*x,2*y-1):
+ fi:
+ end:
+ restriction := array(0..3,0..n_functions-1, 0..n_functions-1):
+ for child from 0 to 3 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][1],
+ support_points[j][2]):
+ od:
+ od:
+ od:
+
+
+ print ("Computing local mass matrix"):
+ # tphi are the basis functions of the linear element. These functions
+ # are used for the computation of the subparametric transformation from
+ # unit cell to real cell.
+ # x and y are arrays holding the x- and y-values of the four vertices
+ # of this cell in real space.
+ x := array(0..3);
+ y := array(0..3);
+ tphi[0] := (1-xi)*(1-eta):
+ tphi[1] := xi*(1-eta):
+ tphi[2] := xi*eta:
+ tphi[3] := (1-xi)*eta:
+ x_real := sum(x[s]*tphi[s], s=0..3):
+ y_real := sum(y[s]*tphi[s], s=0..3):
+ detJ := diff(x_real,xi)*diff(y_real,eta) - diff(x_real,eta)*diff(y_real,xi):
+ for i from 0 to n_functions-1 do
+ print ("line", i):
+ for j from 0 to n_functions-1 do
+ local_mass_matrix[i,j] := int(int(phi_polynom[i] * phi_polynom[j] * detJ,
+ xi=0..1), eta=0..1):
+ od:
+ od:
+
+ print ("computing support points in real space"):
+ for i from 0 to n_functions-1 do
+ real_points[i,0] := subs(xi=support_points[i][1],
+ eta=support_points[i][2], x_real);
+ real_points[i,1] := subs(xi=support_points[i][1],
+ eta=support_points[i][2], y_real);
+ od:
+
+ print ("computing interface constraint matrices"):
+ # compute the interface constraint matrices. these are the values of the
+ # basis functions on the coarse cell's face at the points of the child
+ # cell's basis functions on the child faces
+ face_phi_polynom := array(0..n_face_functions-1):
+ for i from 0 to n_face_functions-1 do
+ # note that the interp function wants vectors indexed from
+ # one and not from zero.
+ values := array(1..n_face_functions):
+ for j from 1 to n_face_functions do
+ values[j] := 0:
+ od:
+ values[i+1] := 1:
+
+ shifted_face_support_points := array (1..n_face_functions):
+ for j from 1 to n_face_functions do
+ shifted_face_support_points[j] := face_support_points[j-1]:
+ od:
+
+ face_phi_polynom[i] := interp (shifted_face_support_points, values, xi):
+ od:
+
+ for i from 0 to 2*(n_face_functions-2)+1-1 do
+ for j from 0 to n_face_functions-1 do
+ interface_constraints[i,j] := subs(xi=constrained_face_support_points[i],
+ face_phi_polynom[j]);
+ od:
+ od:
--- /dev/null
+# --------------------------------- For 2d ---------------------------------
+# -- Use the following maple script to generate the basis functions,
+# -- gradients and prolongation matrices as well as the mass matrix.
+# -- Make sure that the files do not exists beforehand, since output
+# -- is appended instead of overwriting previous contents.
+# --
+# -- You should only have to change the very first lines for polynomials
+# -- of higher order.
+# --------------------------------------------------------------------------
+#
+# $Id$
+# Author: Wolfgang Bangerth, 1998
+
+ n_functions := 16:
+ n_face_functions := 4:
+
+ trial_function := (a1 + a2*xi + a3*xi*xi + a4*xi*xi*xi) +
+ (b1 + b2*xi + b3*xi*xi + b4*xi*xi*xi)*eta +
+ (c1 + c2*xi + c3*xi*xi + c4*xi*xi*xi)*eta*eta +
+ (d1 + d2*xi + d3*xi*xi + d4*xi*xi*xi)*eta*eta*eta:
+ face_trial_function := a + b*xi + c*xi*xi + d*xi*xi*xi:
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ support_points := array(0..n_functions-1):
+ support_points[0] := [0,0]:
+ support_points[1] := [1,0]:
+ support_points[2] := [1,1]:
+ support_points[3] := [0,1]:
+ support_points[4] := [1/3,0]:
+ support_points[5] := [2/3,0]:
+ support_points[6] := [1,1/3]:
+ support_points[7] := [1,2/3]:
+ support_points[8] := [1/3,1]:
+ support_points[9] := [2/3,1]:
+ support_points[10]:= [0,1/3]:
+ support_points[11]:= [0,2/3]:
+ support_points[12]:= [1/3,1/3]:
+ support_points[13]:= [2/3,1/3]:
+ support_points[14]:= [2/3,2/3]:
+ support_points[15]:= [1/3,2/3]:
+
+ face_support_points := array(0..n_face_functions-1):
+ face_support_points[0] := 0:
+ face_support_points[1] := 1:
+ face_support_points[2] := 1/3:
+ face_support_points[3] := 2/3:
+ constrained_face_support_points := array(0..2*(n_face_functions-2)+1-1):
+ constrained_face_support_points[0] := 1/2:
+ constrained_face_support_points[1] := 1/6:
+ constrained_face_support_points[2] := 2/6:
+ constrained_face_support_points[3] := 4/6:
+ constrained_face_support_points[4] := 5/6:
+
+ # do the real work
+ read "lagrange"
+
+
+ # write data to files
+ print ("writing data to files"):
+ readlib(C):
+ C(phi_polynom, filename=cubic2d.shape_value):
+ C(grad_phi_polynom, filename=cubic2d.shape_grad):
+ C(prolongation, filename=cubic2d.prolongation):
+ C(restriction, filename=cubic2d.restriction):
+ C(local_mass_matrix, optimized, filename=cubic2d.massmatrix):
+ C(interface_constraints, filename=cubic2d.constraints):
+ C(real_points, optimized, filename=cubic2d.real_points):
+
--- /dev/null
+# --------------------------------- For 2d ---------------------------------
+# -- Use the following maple script to generate the basis functions,
+# -- gradients and prolongation matrices as well as the mass matrix.
+# -- Make sure that the files do not exists beforehand, since output
+# -- is appended instead of overwriting previous contents.
+# --
+# -- You should only have to change the very first lines for polynomials
+# -- of higher order.
+# --
+# -- Please note:
+# -- Apart from the restriction matrices, I did not initially use it; it is
+# -- an adaption of the script for cubic and quartic elements. For
+# -- some of the data, however, a smaller script is given in the
+# -- FEQuadratic<2> constructor.
+#
+# --------------------------------------------------------------------------
+#
+# $Id$
+# Author: Wolfgang Bangerth, 1998
+
+ n_functions := 9:
+ n_face_functions := 3:
+
+ trial_function := (a1 + a2*xi + a3*xi*xi) +
+ (b1 + b2*xi + b3*xi*xi)*eta +
+ (c1 + c2*xi + c3*xi*xi)*eta*eta:
+ face_trial_function := a + b*xi + c*xi*xi:
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ support_points := array(0..n_functions-1):
+ support_points[0] := [0,0]:
+ support_points[1] := [1,0]:
+ support_points[2] := [1,1]:
+ support_points[3] := [0,1]:
+ support_points[4] := [1/2,0]:
+ support_points[5] := [1,1/2]:
+ support_points[6] := [1/2,1]:
+ support_points[7] := [0,1/2]:
+ support_points[8] := [1/2,1/2]:
+
+ face_support_points := array(0..n_face_functions-1):
+ face_support_points[0] := 0:
+ face_support_points[1] := 1:
+ face_support_points[2] := 1/2:
+
+ constrained_face_support_points := array(0..2*(n_face_functions-2)+1-1):
+ constrained_face_support_points[0] := 1/2:
+ constrained_face_support_points[1] := 1/4:
+ constrained_face_support_points[2] := 3/4:
+
+
+ # do the real work
+ read "lagrange"
+
+
+ # write data to files
+ print ("writing data to files"):
+ readlib(C):
+ C(phi_polynom, filename=quadratic2d.shape_value):
+ C(grad_phi_polynom, filename=quadratic2d.shape_grad):
+ C(prolongation, filename=quadratic2d.prolongation):
+ C(restriction, filename=quadratic2d.restriction):
+ C(local_mass_matrix, optimized, filename=quadratic2d.massmatrix):
+ C(interface_constraints, filename=quadratic2d.constraints):
+ C(real_points, optimized, filename=quadratic2d.real_points):
+
--- /dev/null
+# --------------------------------- For 2d ---------------------------------
+# -- Use the following maple script to generate the basis functions,
+# -- gradients and prolongation matrices as well as the mass matrix.
+# -- Make sure that the files do not exists beforehand, since output
+# -- is appended instead of overwriting previous contents.
+# --
+# -- You should only have to change the very first lines for polynomials
+# -- of higher order.
+# --------------------------------------------------------------------------
+#
+# $Id$
+# Author: Wolfgang Bangerth, 1998
+
+ n_functions := 25:
+ n_face_functions := 5:
+
+ trial_function := (a1 + a2*xi + a3*xi*xi + a4*xi**3 + a5*xi**4) +
+ (b1 + b2*xi + b3*xi*xi + b4*xi**3 + b5*xi**4)*eta +
+ (c1 + c2*xi + c3*xi*xi + c4*xi**3 + c5*xi**4)*eta*eta +
+ (d1 + d2*xi + d3*xi*xi + d4*xi**3 + d5*xi**4)*eta**3 +
+ (e1 + e2*xi + e3*xi*xi + e4*xi**3 + e5*xi**4)*eta**4:
+ face_trial_function := a + b*xi + c*xi*xi + d*xi**3 + e*xi**4:
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ support_points := array(0..n_functions-1):
+ support_points[0] := [0,0]:
+ support_points[1] := [1,0]:
+ support_points[2] := [1,1]:
+ support_points[3] := [0,1]:
+ support_points[4] := [1/4,0]:
+ support_points[5] := [2/4,0]:
+ support_points[6] := [3/4,0]:
+ support_points[7] := [1,1/4]:
+ support_points[8] := [1,2/4]:
+ support_points[9] := [1,3/4]:
+ support_points[10] := [1/4,1]:
+ support_points[11] := [2/4,1]:
+ support_points[12] := [3/4,1]:
+ support_points[13] := [0,1/4]:
+ support_points[14] := [0,2/4]:
+ support_points[15] := [0,3/4]:
+ support_points[16] := [1/4,1/4]:
+ support_points[17] := [3/4,1/4]:
+ support_points[18] := [3/4,3/4]:
+ support_points[19] := [1/4,3/4]:
+ support_points[20] := [1/2,1/4]:
+ support_points[21] := [3/4,1/2]:
+ support_points[22] := [1/2,3/4]:
+ support_points[23] := [1/4,1/2]:
+ support_points[24] := [1/2,1/2]:
+
+ face_support_points := array(0..n_face_functions-1):
+ face_support_points[0] := 0:
+ face_support_points[1] := 1:
+ face_support_points[2] := 1/4:
+ face_support_points[3] := 2/4:
+ face_support_points[4] := 3/4:
+ constrained_face_support_points := array(0..2*(n_face_functions-2)+1-1):
+ constrained_face_support_points[0] := 1/2:
+ constrained_face_support_points[1] := 1/8:
+ constrained_face_support_points[2] := 2/8:
+ constrained_face_support_points[3] := 3/8:
+ constrained_face_support_points[4] := 5/8:
+ constrained_face_support_points[5] := 6/8:
+ constrained_face_support_points[6] := 7/8:
+
+
+ # do the real work
+ read "lagrange"
+
+
+ # write data to files
+ print ("writing data to files"):
+ readlib(C):
+ C(phi_polynom, filename=quartic2d.shape_value):
+ C(grad_phi_polynom, filename=quartic2d.shape_grad):
+ C(prolongation, filename=quartic2d.prolongation):
+ C(restriction, filename=quartic2d.restriction):
+ C(local_mass_matrix, optimized, filename=quartic2d.massmatrix):
+ C(interface_constraints, filename=quartic2d.constraints):
+ C(real_points, optimized, filename=quartic2d.real_points):
+
--- /dev/null
+# Use the following perl scripts to convert the output into a
+# suitable format.
+
+perl -pi -e 's/phi_polynom\[(\d+)\] =/case $1: return/g;' *2d.shape_value
+perl -pi -e 's/([^;])\n/$1/g;' *2d.shape_grad
+perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[0\] = (.*);/case $1: return Point<2>($2,/g;' *2d.shape_grad
+perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[1\] = (.*);/$2);/g;' *2d.shape_grad
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]/($1,$2)/g;' *2d.massmatrix
+perl -pi -e 's/(t\d+) =/const double $1 =/g;' *2d.massmatrix
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' *2d.prolongation
+perl -pi -e 's/.*= 0.0;\n//g;' *2d.prolongation
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' *2d.restriction
+perl -pi -e 's/.*= 0.0;\n//g;' *2d.restriction
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]/($1,$2)/g;' *2d.constraints