--- /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 := 64:
+ n_face_functions := 16:
+ n_constraints := 45:
+
+ 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) +
+ ((e1 + e2*xi + e3*xi*xi + e4*xi*xi*xi) +
+ (f1 + f2*xi + f3*xi*xi + f4*xi*xi*xi)*eta +
+ (g1 + g2*xi + g3*xi*xi + g4*xi*xi*xi)*eta*eta +
+ (h1 + h2*xi + h3*xi*xi + h4*xi*xi*xi)*eta*eta*eta)*zeta +
+ ((i1 + i2*xi + i3*xi*xi + i4*xi*xi*xi) +
+ (j1 + j2*xi + j3*xi*xi + j4*xi*xi*xi)*eta +
+ (k1 + k2*xi + k3*xi*xi + k4*xi*xi*xi)*eta*eta +
+ (l1 + l2*xi + l3*xi*xi + l4*xi*xi*xi)*eta*eta*eta)*zeta*zeta +
+ ((m1 + m2*xi + m3*xi*xi + m4*xi*xi*xi) +
+ (n1 + n2*xi + n3*xi*xi + n4*xi*xi*xi)*eta +
+ (o1 + o2*xi + o3*xi*xi + o4*xi*xi*xi)*eta*eta +
+ (p1 + p2*xi + p3*xi*xi + p4*xi*xi*xi)*eta*eta*eta)*zeta*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, 2, support_points):
+ support_points_fill_quads (32, 2, support_points):
+ support_points_fill_hex (56, 2, support_points):
+
+ 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/3,0]:
+ face_support_points[5] := [2/3,0]:
+ face_support_points[6] := [1,1/3]:
+ face_support_points[7] := [1,2/3]:
+ face_support_points[8] := [1/3,1]:
+ face_support_points[9] := [2/3,1]:
+ face_support_points[10] := [0,1/3]:
+ face_support_points[11] := [0,2/3]:
+ face_support_points[12] := [1/3,1/3]:
+ face_support_points[13] := [1/3,2/3]:
+ face_support_points[14] := [2/3,1/3]:
+ face_support_points[15] := [2/3,2/3]:
+
+ # 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] := 12:
+ constrained_face_function[9] := 13:
+ constrained_face_function[10] := 14:
+ constrained_face_function[11] := 15:
+ constrained_face_function[12] := 32:
+ constrained_face_function[13] := 33:
+ constrained_face_function[14] := 34:
+ constrained_face_function[15] := 35:
+
+ 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/6]): # lines from center to boundary
+ constrained_face_support_points[6] := array(0..1, [1/2,2/6]):
+ constrained_face_support_points[7] := array(0..1, [4/6,1/2]):
+ constrained_face_support_points[8] := array(0..1, [5/6,1/2]):
+ constrained_face_support_points[9] := array(0..1, [1/2,4/6]):
+ constrained_face_support_points[10] := array(0..1, [1/2,5/6]):
+ constrained_face_support_points[11] := array(0..1, [1/6,1/2]):
+ constrained_face_support_points[12] := array(0..1, [2/6,1/2]):
+ constrained_face_support_points[13] := array(0..1, [1/6,0]): # children of bounding lines
+ constrained_face_support_points[14] := array(0..1, [2/6,0]):
+ constrained_face_support_points[15] := array(0..1, [4/6,0]):
+ constrained_face_support_points[16] := array(0..1, [5/6,0]):
+
+ constrained_face_support_points[17] := array(0..1, [1,1/6]):
+ constrained_face_support_points[18] := array(0..1, [1,2/6]):
+ constrained_face_support_points[19] := array(0..1, [1,4/6]):
+ constrained_face_support_points[20] := array(0..1, [1,5/6]):
+
+ constrained_face_support_points[21] := array(0..1, [1/6,1]):
+ constrained_face_support_points[22] := array(0..1, [2/6,1]):
+ constrained_face_support_points[23] := array(0..1, [4/6,1]):
+ constrained_face_support_points[24] := array(0..1, [5/6,1]):
+
+ constrained_face_support_points[25] := array(0..1, [0,1/6]):
+ constrained_face_support_points[26] := array(0..1, [0,2/6]):
+ constrained_face_support_points[27] := array(0..1, [0,4/6]):
+ constrained_face_support_points[28] := array(0..1, [0,5/6]):
+
+ constrained_face_support_points[29] := array(0..1, [1/6,1/6]): # child quads
+ constrained_face_support_points[30] := array(0..1, [2/6,1/6]):
+ constrained_face_support_points[31] := array(0..1, [1/6,2/6]):
+ constrained_face_support_points[32] := array(0..1, [2/6,2/6]):
+
+ constrained_face_support_points[33] := array(0..1, [4/6,1/6]):
+ constrained_face_support_points[34] := array(0..1, [5/6,1/6]):
+ constrained_face_support_points[35] := array(0..1, [4/6,2/6]):
+ constrained_face_support_points[36] := array(0..1, [5/6,2/6]):
+
+ constrained_face_support_points[37] := array(0..1, [4/6,4/6]):
+ constrained_face_support_points[38] := array(0..1, [5/6,4/6]):
+ constrained_face_support_points[39] := array(0..1, [4/6,5/6]):
+ constrained_face_support_points[40] := array(0..1, [5/6,5/6]):
+
+ constrained_face_support_points[41] := array(0..1, [1/6,4/6]):
+ constrained_face_support_points[42] := array(0..1, [2/6,4/6]):
+ constrained_face_support_points[43] := array(0..1, [1/6,5/6]):
+ constrained_face_support_points[44] := array(0..1, [2/6,5/6]):
+
+
+ # do the real work
+ read "lagrange":
+
+
+
+ # write data to files
+ print ("writing data to files"):
+ readlib(C):
+ C(phi_polynom, filename=cubic3d_shape_value):
+ C(grad_phi_polynom, filename=cubic3d_shape_grad):
+ C(grad_grad_phi_polynom, filename=cubic3d_shape_grad_grad):
+ C(prolongation, filename=cubic3d_prolongation):
+ C(restriction, filename=cubic3d_restriction):
+ C(interface_constraints, filename=cubic3d_constraints):
+ C(real_space_points, optimized, filename=cubic3d_real_points):
+
+ writeto (cubic3d_unit_support_points):
+ print (support_points):
+
+
\ No newline at end of file