--- /dev/null
+ dim:=1;
+
+ print (`Computing basis functions`);
+ phi_polynom := array(0..n_functions-1);
+ 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:
+
+ 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);
+ od:
+
+ phi:= proc(i,x,y) 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;
+
+ # find the prolongation matrices such that
+ # phi(k,x,y)|_K_i=prol[i,j,k] child_phi[i](j,x,y)
+ print (`Computing prolongation matrices`):
+ 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;
+
+ # assemble the local mass matrix (on [0,1])
+ # m[i,j]=int_{0..1} phi[i]*phi[j] dx
+ m := array(1..n_functions, 1..n_functions):
+ print (`Assembling mass matrix`):
+ for i from 1 to n_functions do
+ for j from 1 to n_functions do
+ m[i,j] := int(phi_polynom[i-1] * phi_polynom[j-1], xi=0..1);
+ od:
+ od:
+
+ print(`m=`, m);
+
+ # assemble the local mass matrix for child cell 0
+ # m[i,j]=int_{0..0.5}child_phi[0]*child_phi[0] dx
+ child_m := array(1..n_functions, 1..n_functions):
+ child_m:=linalg[scalarmul](m, 1/2**dim);
+
+ print(`Ausgabe=`);
+ print(`child_m=`,child_m);
+
+ # inverte the local mass matrix
+ inv_m := linalg[inverse](m):
+ print(`inv_m=`, inv_m);
+
+ # assembling restriction matrices
+ restriction := array(0..1, 0..n_functions-1, 0..n_functions-1):
+ restr_child := array(1..n_functions, 1..n_functions):
+ prol_child:= array(1..n_functions, 1..n_functions):
+ for child from 0 to 1 do
+ print(`child=`, child);
+ # copy the prologation matrix with a shift 1 and take the transpose
+ for i from 1 to n_functions do
+ for j from 1 to n_functions do
+ prol_child[i,j] := prolongation[child,j-1,i-1]:
+ od:
+ od:
+ restr_child := linalg[multiply](inv_m, prol_child, child_m);
+ print(restr_child);
+ # copy the restriction of this child with a shift 1
+ for i from 1 to n_functions do
+ for j from 1 to n_functions do
+ restriction[child,i-1,j-1] := restr_child[i,j]:
+ od:
+ od:
+ od:
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+# --------------------------------- For 1d ---------------------------------
+# -- Use the following maple script to generate the restriction matrices
+# -- for DG.
+# -- 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: Ralf Hartmann, 2000
+
+# for DG(1)
+
+ n_functions := 2:
+
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ # phi(i,support_points[j])=delta_ij
+ support_points := array(0..n_functions-1):
+ support_points[0] := 0:
+ support_points[1] := 1:
+
+ read restriction_dg;
+
+ print (`writing data to files`):
+ readlib(C):
+ C(restriction, filename=restriction_dg1_txt);
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+# --------------------------------- For 1d ---------------------------------
+# -- Use the following maple script to generate the restriction matrices
+# -- for DG.
+# -- 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: Ralf Hartmann, 2000
+
+# for DG(2)
+
+ n_functions := 3:
+
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ # phi(i,support_points[j])=delta_ij
+ support_points := array(0..n_functions-1):
+ support_points[0] := 0:
+ support_points[1] := 1:
+ support_points[2] := 1/2:
+
+ read restriction_dg;
+
+ print (`writing data to files`):
+ readlib(C):
+ C(restriction, filename=restriction_dg2_txt);
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+# --------------------------------- For 1d ---------------------------------
+# -- Use the following maple script to generate the restriction matrices
+# -- for DG.
+# -- 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: Ralf Hartmann, 2000
+
+# for DG(3)
+
+ n_functions := 4:
+
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ # phi(i,support_points[j])=delta_ij
+ 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;
+
+ read restriction_dg;
+
+ print (`writing data to files`):
+ readlib(C):
+ C(restriction, filename=restriction_dg3_txt);
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+# --------------------------------- For 1d ---------------------------------
+# -- Use the following maple script to generate the restriction matrices
+# -- for DG.
+# -- 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: Ralf Hartmann, 2000
+
+# for DG(4)
+
+ n_functions := 5:
+
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ # phi(i,support_points[j])=delta_ij
+ 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;
+
+ read restriction_dg;
+
+ print (`writing data to files`):
+ readlib(C):
+ C(restriction, filename=restriction_dg4_txt);
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+ dim:=2;
+
+ print (`Computing basis functions`);
+ phi_polynom := array(0..n_functions-1);
+ 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);
+ 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
+ # child_phi[c](i, points[c][j, ])=delta_ij
+ 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:
+
+ # find the prolongation matrices such that
+ # phi(k,x,y)|_K_i=prol[i,j,k] child_phi[i](j,x,y)
+ print (`Computing prolongation matrices`):
+ prolongation := array(0..3,0..n_functions-1, 0..n_functions-1):
+ for i from 0 to 3 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,1], points[i][j,2]);
+ od:
+ od:
+ od:
+
+ # assemble the local mass matrix (on the unit square)
+ # m[i,j]=int_{0..1}int_{0..1} phi[i]*phi[j] dxdy
+ m := array(1..n_functions, 1..n_functions):
+ print (`Assembling mass matrix`):
+ for i from 1 to n_functions do
+ for j from 1 to n_functions do
+ m[i,j] := int(int(phi_polynom[i-1] * phi_polynom[j-1], xi=0..1), eta=0..1);
+ od:
+ od:
+
+ print(`m=`, m);
+
+ # assemble the local mass matrix for child cell 0
+ # m[i,j]=int_{0..0.5}int_{0..0.5} child_phi[0]*child_phi[0] dxdy
+ child_m := array(1..n_functions, 1..n_functions):
+ child_m:=linalg[scalarmul](m, 1/2**dim);
+
+ print(`Ausgabe=`);
+ print(`child_m=`,child_m);
+
+ # inverte the local mass matrix
+ inv_m := linalg[inverse](m):
+ print(`inv_m=`, inv_m);
+
+ # assembling restriction matrices
+ restriction := array(0..3, 0..n_functions-1, 0..n_functions-1):
+ restr_child := array(1..n_functions, 1..n_functions):
+ prol_child:= array(1..n_functions, 1..n_functions):
+ for child from 0 to 3 do
+ print(`child=`, child);
+ # copy the prologation matrix with a shift 1 and take the transponent
+ for i from 1 to n_functions do
+ for j from 1 to n_functions do
+ prol_child[i,j] := prolongation[child,j-1,i-1]:
+ od:
+ od:
+ restr_child := linalg[multiply](inv_m, prol_child, child_m);
+ print(restr_child);
+ # copy the restriction of this child with a shift 1
+ for i from 1 to n_functions do
+ for j from 1 to n_functions do
+ restriction[child,i-1,j-1] := restr_child[i,j]:
+ od:
+ od:
+ od:
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+# --------------------------------- For 2d ---------------------------------
+# -- Use the following maple script to generate the restriction matrices
+# -- for DG.
+# -- 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: Ralf Hartmann, 2000
+
+# for DG(1)
+
+ n_functions := 4:
+
+ trial_function := (a1 + a2*xi) +
+ (b1 + b2*xi)*eta:
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ # phi(i,support_points[j])=delta_ij
+ 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]:
+
+ read restriction_dg;
+
+ print (`writing data to files`):
+ readlib(C):
+ C(restriction, filename=restriction_dg1_txt);
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+# --------------------------------- For 2d ---------------------------------
+# -- Use the following maple script to generate the restriction matrices
+# -- for DG.
+# -- 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: Ralf Hartmann, 2000
+
+# for DG(2)
+
+ n_functions := 9:
+
+ trial_function := (a1 + a2*xi + a3*xi*xi) +
+ (b1 + b2*xi + b3*xi*xi)*eta +
+ (c1 + c2*xi + c3*xi*xi)*eta*eta:
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ # phi(i,support_points[j])=delta_ij
+ 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]:
+
+ read restriction_dg;
+
+ print (`writing data to files`):
+ readlib(C):
+ C(restriction, filename=restriction_dg2_txt);
+
--- /dev/null
+# --------------------------------- For 2d ---------------------------------
+# -- Use the following maple script to generate the restriction matrices
+# -- for DG.
+# -- 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: Ralf Hartmann, 2000
+
+# for DG(3)
+
+ n_functions := 16:
+
+ 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:
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ # phi(i,support_points[j])=delta_ij
+ 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]:
+
+ read restriction_dg;
+
+ print (`writing data to files`):
+ readlib(C):
+ C(restriction, filename=restriction_dg3_txt);
+
+
+
+
+
+
+
+
+
+
+
--- /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: Ralf Hartmann, 2000
+
+# for DG(4)
+
+ 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:
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ # phi(i,support_points[j])=delta_ij
+ 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]:
+
+ read restriction_dg;
+
+ print (`writing data to files`):
+ readlib(C):
+ C(restriction, filename=restriction_dg4_txt);
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+ dim:=3;
+
+ print (`Computing basis functions`);
+ phi_polynom := array(0..n_functions-1);
+ 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);
+ 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
+ # child_phi[c](i, points[c][j, ])=delta_ij
+ 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:
+
+ # find the prolongation matrices such that
+ # phi(k,x,y,z)|_K_i=prol[i,j,k] child_phi[i](j,x,y,z)
+ print (`Computing prolongation matrices`):
+ prolongation := array(0..7,0..n_functions-1, 0..n_functions-1):
+ for i from 0 to 7 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,1], points[i][j,2], points[i][j,3]);
+ od:
+ od:
+ od:
+
+ # assemble the local mass matrix (on the unit square)
+ # m[i,j]=int_{0..1}int_{0..1}int_{0..1} phi[i]*phi[j] dxdydz
+ m := array(1..n_functions, 1..n_functions):
+ print (`Assembling mass matrix`):
+ for i from 1 to n_functions do
+ for j from 1 to n_functions do
+ m[i,j] := int(int(int(phi_polynom[i-1] * phi_polynom[j-1], xi=0..1), eta=0..1), zeta=0..1);
+ od:
+ od:
+
+ print(`m=`, m);
+
+ # assemble the local mass matrix for child cell 0
+ # m[i,j]=int_{0..0.5}int_{0..0.5}int_{0..0.5} child_phi[0]*child_phi[0] dxdydz
+ child_m := array(1..n_functions, 1..n_functions):
+ child_m:=linalg[scalarmul](m, 1/2**dim);
+
+ print(`Ausgabe=`);
+ print(`child_m=`,child_m);
+
+ # inverte the local mass matrix
+ inv_m := linalg[inverse](m):
+ print(`inv_m=`, inv_m);
+
+ # assembling restriction matrices
+ restriction := array(0..7, 0..n_functions-1, 0..n_functions-1):
+ restr_child := array(1..n_functions, 1..n_functions):
+ prol_child:= array(1..n_functions, 1..n_functions):
+ for child from 0 to 7 do
+ print(`child=`, child);
+ # copy the prologation matrix with a shift 1 and take the transponent
+ for i from 1 to n_functions do
+ for j from 1 to n_functions do
+ prol_child[i,j] := prolongation[child,j-1,i-1]:
+ od:
+ od:
+ restr_child := linalg[multiply](inv_m, prol_child, child_m);
+ print(restr_child);
+ # copy the restriction of this child with a shift 1
+ for i from 1 to n_functions do
+ for j from 1 to n_functions do
+ restriction[child,i-1,j-1] := restr_child[i,j]:
+ od:
+ od:
+ od:
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+# --------------------------------- For 3d ---------------------------------
+# -- Use the following maple script to generate the restriction matrices
+# -- for DG.
+# -- 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: Ralf Hartmann, 2000
+
+# for DG(1)
+
+ n_functions := 8:
+
+ trial_function := ((a1 + a2*xi) +
+ (b1 + b2*xi)*eta) +
+ ((d1 + d2*xi) +
+ (e1 + e2*xi)*eta)*zeta:
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ # phi(i,support_points[j])=delta_ij
+ 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]):
+
+ read restriction_dg;
+
+ print (`writing data to files`):
+ readlib(C):
+ C(restriction, filename=restriction_dg1_txt);
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+# --------------------------------- For 3d ---------------------------------
+# -- Use the following maple script to generate the restriction matrices
+# -- for DG.
+# -- 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: Ralf Hartmann, 2000
+
+# for DG(2)
+
+ read lagrange_tools:
+
+ n_functions := 27:
+
+ 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:
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ # phi(i,support_points[j])=delta_ij
+ 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
+
+ read restriction_dg;
+
+ print (`writing data to files`):
+ readlib(C):
+ C(restriction, filename=restriction_dg2_txt);
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+# --------------------------------- For 3d ---------------------------------
+# -- Use the following maple script to generate the restriction matrices
+# -- for DG.
+# -- 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: Ralf Hartmann, 2000
+
+# for DG(3)
+
+ read lagrange_tools:
+
+ n_functions := 64:
+
+ 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:
+ # note: support_points[i] is a vector which is indexed from
+ # one and not from zero!
+ # phi(i,support_points[j])=delta_ij
+ 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):
+
+ read restriction_dg;
+
+ print (`writing data to files`):
+ readlib(C):
+ C(restriction, filename=restriction_dg3_txt);
+
+
+
+
+
+
+
+
+
+
+
projects / dealii-svn.git / commitdiff