--- /dev/null
+ # 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);
+ z := 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):
+ z_real := sum(z[s]*tphi[s], s=0..3):
+
+ image := vector([x_real, y_real, z_real]):
+
+ outward_vector := linalg[crossprod] (map(diff, image, xi),
+ map(diff, image,eta)):
+ detJ := linalg[norm] (outward_vector, 2):
+ normal_vector := map (p->p/detJ, outward_vector):
+
+
+ readlib(C):
+ C(detJ, optimized):
+
+ # apply the following perl scripts:
+ # perl -pi -e 's/x\[(\d)\]/vertices[$1](0)/g; s/y\[(\d)\]/vertices[$1](1)/g; s/z\[(\d)\]/vertices[$1](2)/g;'
+ # perl -pi -e 's/^\s*t/const double t/g;'
\ No newline at end of file
--- /dev/null
+# Maple script to compute much of the data needed to implement the
+# family of Lagrange elements in 3d. 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,0..2):
+ grad_grad_phi_polynom := array(0..n_functions-1,0..2,0..2):
+ local_mass_matrix := array(0..n_functions-1, 0..n_functions-1):
+ prolongation := array(0..7,0..n_functions-1, 0..n_functions-1):
+ interface_constraints := array(0..n_constraints-1,
+ 0..n_face_functions-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],
+ zeta=support_points[j][3],
+ 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):
+ grad_phi_polynom[i,2] := diff(phi_polynom[i], zeta):
+
+ grad_grad_phi_polynom[i,0,0] := diff(phi_polynom[i], xi, xi):
+ grad_grad_phi_polynom[i,0,1] := diff(phi_polynom[i], xi, eta):
+ grad_grad_phi_polynom[i,0,2] := diff(phi_polynom[i], xi, zeta):
+ grad_grad_phi_polynom[i,1,0] := diff(phi_polynom[i], eta,xi):
+ grad_grad_phi_polynom[i,1,1] := diff(phi_polynom[i], eta,eta):
+ grad_grad_phi_polynom[i,1,2] := diff(phi_polynom[i], eta,zeta):
+ grad_grad_phi_polynom[i,2,0] := diff(phi_polynom[i], zeta,xi):
+ grad_grad_phi_polynom[i,2,1] := diff(phi_polynom[i], zeta,eta):
+ grad_grad_phi_polynom[i,2,2] := diff(phi_polynom[i], zeta,zeta):
+ od:
+
+ phi:= proc(i,x,y,z) subs(xi=x, eta=y, zeta=z, 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..3):
+ points[1] := array(0..n_functions-1, 1..3):
+ points[2] := array(0..n_functions-1, 1..3):
+ points[3] := array(0..n_functions-1, 1..3):
+ points[4] := array(0..n_functions-1, 1..3):
+ points[5] := array(0..n_functions-1, 1..3):
+ points[6] := array(0..n_functions-1, 1..3):
+ points[7] := array(0..n_functions-1, 1..3):
+ 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[0][i,3] := support_points[i][3]/2:
+
+ points[1][i,1] := support_points[i][1]/2+1/2:
+ points[1][i,2] := support_points[i][2]/2:
+ points[1][i,3] := support_points[i][3]/2:
+
+ points[2][i,1] := support_points[i][1]/2+1/2:
+ points[2][i,2] := support_points[i][2]/2:
+ points[2][i,3] := support_points[i][3]/2+1/2:
+
+ points[3][i,1] := support_points[i][1]/2:
+ points[3][i,2] := support_points[i][2]/2:
+ points[3][i,3] := support_points[i][3]/2+1/2:
+
+ points[4][i,1] := support_points[i][1]/2:
+ points[4][i,2] := support_points[i][2]/2+1/2:
+ points[4][i,3] := support_points[i][3]/2:
+
+ points[5][i,1] := support_points[i][1]/2+1/2:
+ points[5][i,2] := support_points[i][2]/2+1/2:
+ points[5][i,3] := support_points[i][3]/2:
+
+ points[6][i,1] := support_points[i][1]/2+1/2:
+ points[6][i,2] := support_points[i][2]/2+1/2:
+ points[6][i,3] := support_points[i][3]/2+1/2:
+
+ points[7][i,1] := support_points[i][1]/2:
+ points[7][i,2] := support_points[i][2]/2+1/2:
+ points[7][i,3] := support_points[i][3]/2+1/2:
+ od:
+
+ print ("Computing prolongation matrices"):
+ for i from 0 to 7 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], points[i][j,3]):
+ 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, z)
+ if ((x>1/2) or (y>1/2) or (z>1/2)) then
+ 0:
+ else
+ phi(i,2*x,2*y,2*z):
+ fi:
+ end:
+ child_phi[1] := proc(i, x, y, z)
+ if ((x<1/2) or (y>1/2) or (z>1/2)) then
+ 0:
+ else
+ phi(i,2*x-1,2*y, 2*z):
+ fi:
+ end:
+ child_phi[2] := proc(i, x, y, z)
+ if ((x<1/2) or (y>1/2) or (z<1/2)) then
+ 0:
+ else
+ phi(i,2*x-1,2*y, 2*z-1):
+ fi:
+ end:
+ child_phi[3] := proc(i, x, y, z)
+ if ((x>1/2) or (y>1/2) or (z<1/2)) then
+ 0:
+ else
+ phi(i,2*x,2*y,2*z-1):
+ fi:
+ end:
+ child_phi[4] := proc(i, x, y, z)
+ if ((x>1/2) or (y<1/2) or (z>1/2)) then
+ 0:
+ else
+ phi(i,2*x,2*y-1,2*z):
+ fi:
+ end:
+ child_phi[5] := proc(i, x, y, z)
+ if ((x<1/2) or (y<1/2) or (z>1/2)) then
+ 0:
+ else
+ phi(i,2*x-1,2*y-1,2*z):
+ fi:
+ end:
+ child_phi[6] := proc(i, x, y, z)
+ if ((x<1/2) or (y<1/2) or (z<1/2)) then
+ 0:
+ else
+ phi(i,2*x-1,2*y-1,2*z-1):
+ fi:
+ end:
+ child_phi[7] := proc(i, x, y, z)
+ if ((x>1/2) or (y<1/2) or (z<1/2)) then
+ 0:
+ else
+ phi(i,2*x,2*y-1,2*z-1):
+ fi:
+ end:
+ restriction := array(0..7,0..n_functions-1, 0..n_functions-1):
+ for child from 0 to 7 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],
+ support_points[j][3]):
+ od:
+ od:
+ 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
+ #
+ # first compute for each function on the (large) face the polynom
+ # we get this by evaluating the respective global trial function
+ # with y=0
+ face_phi_polynom := array(0..n_face_functions-1):
+ for j from 0 to n_face_functions-1 do
+ face_phi_polynom[j] := proc(xi,eta)
+ subs(dummy=0, phi(constrained_face_function[j],xi,dummy,eta)):
+ end:
+ od:
+
+ for i from 0 to n_constraints-1 do
+ for j from 0 to n_face_functions-1 do
+ interface_constraints[i,j]
+ := face_phi_polynom[j](constrained_face_support_points[i][0],
+ constrained_face_support_points[i][1]):
+ od:
+ od:
+
+
+ # 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. same for z
+ #
+ print ("Computing real space support points"):
+ x := array(0..7);
+ y := array(0..7);
+ z := array(0..7):
+ tphi[0] := (1-xi)*(1-eta)*(1-zeta):
+ tphi[1] := xi*(1-eta)*(1-zeta):
+ tphi[2] := xi*(1-eta)*zeta:
+ tphi[3] := (1-xi)*(1-eta)*zeta:
+ tphi[4] := (1-xi)*eta*(1-zeta):
+ tphi[5] := xi*eta*(1-zeta):
+ tphi[6] := xi*eta*zeta:
+ tphi[7] := (1-xi)*eta*zeta:
+ x_real := sum(x[s]*tphi[s], s=0..7):
+ y_real := sum(y[s]*tphi[s], s=0..7):
+ z_real := sum(z[s]*tphi[s], s=0..7):
+
+ real_space_points := array(0..n_functions-1,0..2):
+ for i from 0 to n_functions-1 do
+ real_space_points[i,0] :=
+ subs(xi=support_points[i][1],
+ eta=support_points[i][2],
+ zeta=support_points[i][3],
+ x_real):
+ real_space_points[i,1] :=
+ subs(xi=support_points[i][1],
+ eta=support_points[i][2],
+ zeta=support_points[i][3],
+ y_real):
+ real_space_points[i,2] :=
+ subs(xi=support_points[i][1],
+ eta=support_points[i][2],
+ zeta=support_points[i][3],
+ z_real):
+ od:
\ No newline at end of file
--- /dev/null
+# --------------------------------- For 3d ---------------------------------
+# -- 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, 1999
+
+ n_functions := 8:
+ n_face_functions := 4:
+ n_constraints := 5:
+
+ trial_function := ((a1 + a2*xi) +
+ (b1 + b2*xi)*eta) +
+ ((d1 + d2*xi) +
+ (e1 + e2*xi)*eta)*zeta:
+ face_trial_function := subs(zeta=0, trial_function):
+ # 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] := array(1..3, [0,0,0]):
+ support_points[1] := array(1..3, [1,0,0]):
+ support_points[2] := array(1..3, [1,0,1]):
+ support_points[3] := array(1..3, [0,0,1]):
+ support_points[4] := array(1..3, [0,1,0]):
+ support_points[5] := array(1..3, [1,1,0]):
+ support_points[6] := array(1..3, [1,1,1]):
+ support_points[7] := array(1..3, [0,1,1]):
+
+ face_support_points := array(0..n_face_functions-1):
+ face_support_points[0] := [0,0]:
+ face_support_points[1] := [1,0]:
+ face_support_points[2] := [1,1]:
+ face_support_points[3] := [0,1]:
+
+ # list of functions which are at face 0, used to compute
+ # the constraints on a face
+ constrained_face_function := array (0..n_face_functions-1):
+ # the list of points at which we want the functions at
+ # faces to be evaluated
+ constrained_face_support_points := array(0..n_constraints-1):
+ constrained_face_function[0] := 0:
+ constrained_face_function[1] := 1:
+ constrained_face_function[2] := 2:
+ constrained_face_function[3] := 3:
+ constrained_face_support_points[0] := array(0..1, [1/2,1/2]):
+ constrained_face_support_points[1] := array(0..1, [1/2,0]):
+ constrained_face_support_points[2] := array(0..1, [1,1/2]):
+ constrained_face_support_points[3] := array(0..1, [1/2,1]):
+ constrained_face_support_points[4] := array(0..1, [0,1/2]):
+
+
+ # do the real work
+ read "lagrange":
+
+ # ... originally taken from another comment, so this does not
+ # fit in here too well...
+ #
+ # Since we're already at it and need it anyway, we also compute the
+ # Jacobian matrix of the transform and its derivatives. For the
+ # question of whether to take the given form or its transpose, refer
+ # to the documentation of the FEValues class and the source code
+ # documentation of FELinearMapping::fill_fe_values. Also note, that
+ # the computed inverse is multiplied to the unit cell gradients
+ # *from the right*.
+ print ("Computing Jacobian matrices"):
+ Jacobian := linalg[matrix](3,3, [[diff(x_real,xi), diff(x_real,eta), diff(x_real,zeta)],
+ [diff(y_real,xi), diff(y_real,eta), diff(y_real,zeta)],
+ [diff(z_real,xi), diff(z_real,eta), diff(z_real,zeta)]]):
+ inverseJacobian := linalg[inverse](Jacobian):
+ detJ := linalg[det](Jacobian):
+
+ grad_inverseJacobian := array(1..3, 1..3, 1..3):
+ for i from 1 to 3 do
+ for j from 1 to 3 do
+ for k from 1 to 3 do
+ if (i=1) then
+ grad_inverseJacobian[i,j,k] := diff (inverseJacobian[j,k], xi):
+ else
+ if (i=2) then
+ grad_inverseJacobian[i,j,k] := diff (inverseJacobian[j,k], eta):
+ else
+ grad_inverseJacobian[i,j,k] := diff (inverseJacobian[j,k], zeta):
+ fi:
+ fi:
+ od:
+ od:
+ od:
+
+
+ print ("computing support points in real space"):
+ real_points := array(0..n_functions-1, 0..2);
+ for i from 0 to n_functions-1 do
+ real_points[i,0] := subs(xi=support_points[i][1],
+ eta=support_points[i][2],
+ zeta=support_points[i][3], x_real);
+ real_points[i,1] := subs(xi=support_points[i][1],
+ eta=support_points[i][2],
+ zeta=support_points[i][3], y_real);
+ real_points[i,2] := subs(xi=support_points[i][1],
+ eta=support_points[i][2],
+ zeta=support_points[i][3], z_real);
+ od:
+
+
+
+ # write data to files
+ print ("writing data to files"):
+ readlib(C):
+ C(phi_polynom, filename=linear3d_shape_value):
+ C(grad_phi_polynom, filename=linear3d_shape_grad):
+ C(grad_grad_phi_polynom, filename=linear3d_shape_grad_grad):
+ C(prolongation, filename=linear3d_prolongation):
+ C(restriction, filename=linear3d_restriction):
+ C(local_mass_matrix, filename=linear3d_massmatrix):
+ C(interface_constraints, filename=linear3d_constraints):
+ C(real_points, optimized, filename=linear3d_real_points):
+ C(inverseJacobian, filename=linear3d_inverse_jacobian):
+ C(grad_inverseJacobian, filename=linear3d_inverse_jacobian_grad):
--- /dev/null
+# --------------------------------- For 3d ---------------------------------
+# -- 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, 1999
+
+ read "lagrange-tools":
+
+ n_functions := 27:
+ n_face_functions := 9:
+ n_constraints := 21:
+
+ trial_function := ((a1 + a2*xi + a3*xi*xi) +
+ (b1 + b2*xi + b3*xi*xi)*eta +
+ (c1 + c2*xi + c3*xi*xi)*eta*eta) +
+ ((d1 + d2*xi + d3*xi*xi) +
+ (e1 + e2*xi + e3*xi*xi)*eta +
+ (f1 + f2*xi + f3*xi*xi)*eta*eta)*zeta +
+ ((g1 + g2*xi + g3*xi*xi) +
+ (h1 + h2*xi + h3*xi*xi)*eta +
+ (i1 + i2*xi + i3*xi*xi)*eta*eta)*zeta*zeta:
+ face_trial_function := subs(zeta=0, trial_function):
+ # 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_fill_vertices (0, support_points):
+ support_points_fill_lines (8, 1, support_points):
+ support_points[20] := array(1..3, [1/2, 0, 1/2]): #faces
+ support_points[21] := array(1..3, [1/2, 1, 1/2]):
+ support_points[22] := array(1..3, [1/2, 1/2, 0]):
+ support_points[23] := array(1..3, [1, 1/2, 1/2]):
+ support_points[24] := array(1..3, [1/2, 1/2, 1]):
+ support_points[25] := array(1..3, [0, 1/2, 1/2]):
+ support_points[26] := array(1..3, [1/2, 1/2,1/2]): #center
+
+ face_support_points := array(0..n_face_functions-1):
+ face_support_points[0] := [0,0]:
+ face_support_points[1] := [1,0]:
+ face_support_points[2] := [1,1]:
+ face_support_points[3] := [0,1]:
+ face_support_points[4] := [1/2,0]:
+ face_support_points[5] := [1,1/2]:
+ face_support_points[6] := [1/2,1]:
+ face_support_points[7] := [0,1/2]:
+ face_support_points[8] := [1/2,1/2]:
+
+ # list of functions which are at face 0, used to compute
+ # the constraints on a face
+ constrained_face_function := array (0..n_face_functions-1):
+ # the list of points at which we want the functions at
+ # faces to be evaluated
+ constrained_face_support_points := array(0..n_constraints-1):
+ constrained_face_function[0] := 0:
+ constrained_face_function[1] := 1:
+ constrained_face_function[2] := 2:
+ constrained_face_function[3] := 3:
+ constrained_face_function[4] := 8:
+ constrained_face_function[5] := 9:
+ constrained_face_function[6] := 10:
+ constrained_face_function[7] := 11:
+ constrained_face_function[8] := 20:
+ constrained_face_support_points[0] := array(0..1, [1/2,1/2]): # center vertex
+ constrained_face_support_points[1] := array(0..1, [1/2,0]): # centers of large lines
+ constrained_face_support_points[2] := array(0..1, [1,1/2]):
+ constrained_face_support_points[3] := array(0..1, [1/2,1]):
+ constrained_face_support_points[4] := array(0..1, [0,1/2]):
+ constrained_face_support_points[5] := array(0..1, [1/2,1/4]): # lines from center to boundary
+ constrained_face_support_points[6] := array(0..1, [3/4,1/2]):
+ constrained_face_support_points[7] := array(0..1, [1/2,3/4]):
+ constrained_face_support_points[8] := array(0..1, [1/4,1/2]):
+ constrained_face_support_points[9] := array(0..1, [1/4,0]): # children of bounding lines
+ constrained_face_support_points[10] := array(0..1, [3/4,0]):
+ constrained_face_support_points[11] := array(0..1, [1,1/4]):
+ constrained_face_support_points[12] := array(0..1, [1,3/4]):
+ constrained_face_support_points[13] := array(0..1, [1/4,1]):
+ constrained_face_support_points[14] := array(0..1, [3/4,1]):
+ constrained_face_support_points[15] := array(0..1, [0,1/4]):
+ constrained_face_support_points[16] := array(0..1, [0,3/4]):
+ constrained_face_support_points[17] := array(0..1, [1/4,1/4]): # child quads
+ constrained_face_support_points[18] := array(0..1, [3/4,1/4]):
+ constrained_face_support_points[19] := array(0..1, [3/4,3/4]):
+ constrained_face_support_points[20] := array(0..1, [1/4,3/4]):
+
+
+ # do the real work
+ read "lagrange":
+
+
+
+ # write data to files
+ print ("writing data to files"):
+ readlib(C):
+ C(phi_polynom, filename=quadratic3d_shape_value):
+ C(grad_phi_polynom, filename=quadratic3d_shape_grad):
+ C(grad_grad_phi_polynom, filename=quadratic3d_shape_grad_grad):
+ C(prolongation, filename=quadratic3d_prolongation):
+ C(restriction, filename=quadratic3d_restriction):
+ C(interface_constraints, filename=quadratic3d_constraints):
+ C(real_space_points, optimized, filename=quadratic3d_real_points):
+
+ writeto (quadratic3d_unit_support_points):
+ print (support_points):
+
+
\ No newline at end of file
--- /dev/null
+support_points_fill_vertices := proc (starting_index, support_points)
+ support_points[starting_index+0] := array(1..3, [0,0,0]):
+ support_points[starting_index+1] := array(1..3, [1,0,0]):
+ support_points[starting_index+2] := array(1..3, [1,0,1]):
+ support_points[starting_index+3] := array(1..3, [0,0,1]):
+ support_points[starting_index+4] := array(1..3, [0,1,0]):
+ support_points[starting_index+5] := array(1..3, [1,1,0]):
+ support_points[starting_index+6] := array(1..3, [1,1,1]):
+ support_points[starting_index+7] := array(1..3, [0,1,1]):
+end:
+
+
+
+support_points_fill_lines := proc (starting_index, dofs_per_line, support_points)
+
+ local next_index, increment, i:
+
+ next_index := starting_index:
+ increment := 1/(dofs_per_line+1):
+
+ # line 0
+ for i from 1 to dofs_per_line do
+ print (next_index+i-1):
+ support_points[next_index+i-1]
+ := array (1..3, [i*increment, 0, 0]):
+ next_index := next_index+1
+ od:
+
+ # line 1
+ for i from 1 to dofs_per_line do
+ support_points[next_index+i-1]
+ := array (1..3, [1, 0, i*increment]):
+ next_index := next_index+1
+ od:
+
+ # line 2
+ for i from 1 to dofs_per_line do
+ support_points[next_index+i-1]
+ := array (1..3, [i*increment, 0, 1]):
+ next_index := next_index+1
+ od:
+
+ # line 3
+ for i from 1 to dofs_per_line do
+ support_points[next_index+i-1]
+ := array (1..3, [0, 0, i*increment]):
+ next_index := next_index+1
+ od:
+
+ # line 4
+ for i from 1 to dofs_per_line do
+ support_points[next_index+i-1]
+ := array (1..3, [i*increment, 1, 0]):
+ next_index := next_index+1
+ od:
+
+ # line 5
+ for i from 1 to dofs_per_line do
+ support_points[next_index+i-1]
+ := array (1..3, [1, 1, i*increment]):
+ next_index := next_index+1
+ od:
+
+ # line 6
+ for i from 1 to dofs_per_line do
+ support_points[next_index+i-1]
+ := array (1..3, [i*increment, 1, 1]):
+ next_index := next_index+1
+ od:
+
+ # line 7
+ for i from 1 to dofs_per_line do
+ support_points[next_index+i-1]
+ := array (1..3, [0, 1, i*increment]):
+ next_index := next_index+1
+ od:
+
+
+ # line 8
+ for i from 1 to dofs_per_line do
+ support_points[next_index+i-1]
+ := array (1..3, [0, i*increment,0]):
+ next_index := next_index+1
+ od:
+
+ # line 9
+ for i from 1 to dofs_per_line do
+ support_points[next_index+i-1]
+ := array (1..3, [1, i*increment, 0]):
+ next_index := next_index+1
+ od:
+
+ # line 10
+ for i from 1 to dofs_per_line do
+ support_points[next_index+i-1]
+ := array (1..3, [1, i*increment, 1]):
+ next_index := next_index+1
+ od:
+
+ # line 11
+ for i from 1 to dofs_per_line do
+ support_points[next_index+i-1]
+ := array (1..3, [0, i*increment, 1]):
+ next_index := next_index+1
+ od:
+end:
\ No newline at end of file
--- /dev/null
+# Use the following perl scripts to convert the output into a
+# suitable format.
+#
+# $Id$
+# Wolfgang Bangerth, 1998
+
+perl -pi -e 's/phi_polynom\[(\d+)\] =/case $1: return/g;' *3d_shape_value
+
+perl -pi -e 's/([^;])\n/$1/g;' *3d_shape_grad
+perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[0\] = (.*);/case $1: return Point<3>($2,/g;' *3d_shape_grad
+perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[[01]\] = (.*);/$2,/g;' *3d_shape_grad
+perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[2\] = (.*);/$2);/g;' *3d_shape_grad
+
+
+# concatenate all lines for each entry
+perl -pi -e 's/([^;])\n/$1/g;' *3d_shape_grad_grad
+# rename the variable
+perl -pi -e 's/\s*grad_grad_phi_polynom/return_value/g;' *3d_shape_grad_grad
+# insert 'case' and 'break' statements
+perl -pi -e 's/(return_value\[(\d)\]\[0\]\[0\] = .*;)/break;\ncase $2:\n$1/g;' *3d_shape_grad_grad
+# eliminate first index, since that one is caught by the 'case' statement
+perl -pi -e 's/return_value\[\d+\]/return_value/g;' *3d_shape_grad_grad
+# delete lines where only a zero is set, since this already is done in the constructor
+perl -pi -e 's/.*= 0.0;\n//g;' *3d_shape_grad_grad
+
+
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' *3d_prolongation
+perl -pi -e 's/.*= 0.0;\n//g;' *3d_prolongation
+
+
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' *3d_restriction
+perl -pi -e 's/.*= 0.0;\n//g;' *3d_restriction
+
+
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]/($1,$2)/g;' *3d_constraints
+perl -pi -e 's/.*= 0.0;\n//g;' *3d_constraints
+
+perl -pi -e 's/^\s*t/const double t/g;' *3d_inverse_jacobian
+perl -pi -e 's/x\[(\d)\]/vertices[$1](0)/g;' *3d_inverse_jacobian
+perl -pi -e 's/y\[(\d)\]/vertices[$1](1)/g;' *3d_inverse_jacobian
+perl -pi -e 's/inverseJacobian/jacobians[point]/g;' *3d_inverse_jacobian
+perl -pi -e 's/\[(\d)\]\[(\d)\] =/($1,$2) =/g;' *3d_inverse_jacobian
+
+perl -pi -e 's/^\s*t/const double t/g;' *3d_inverse_jacobian_grad
+perl -pi -e 's/x\[(\d)\]/vertices[$1](0)/g;' *3d_inverse_jacobian_grad
+perl -pi -e 's/y\[(\d)\]/vertices[$1](1)/g;' *3d_inverse_jacobian_grad
+perl -pi -e 's/inverseJacobian/jacobians_grad[point]/g;' *3d_inverse_jacobian_grad
+
+
+perl -pi -e 's/^array.*\n//g; s/^\s*\]\)//g; s/^\n//g;' *3d_unit_support_points
+perl -pi -e 's/\s+\((\d+)\)/ unit_points[$1]/g;' *3d_unit_support_points
+perl -pi -e 's/= \[/= Point<3>(/g; s/\]\s*\n/);\n/g;' *3d_unit_support_points
+
+
+perl -pi -e 's/real_space_points\[(\d+)\]\[(\d+)\]/support_points[$1]($2)/g;' *3d_real_points
+perl -pi -e 's/x\[(\d+)\]/vertices[$1](0)/g; s/y\[(\d+)\]/vertices[$1](1)/g; s/z\[(\d+)\]/vertices[$1](2)/g;' *3d_real_points
\ No newline at end of file