From: Ralf Hartmann Date: Mon, 22 May 2000 12:38:30 +0000 (+0000) Subject: scripts for assembling the restriction matrices for DG elements. Restriction for... X-Git-Tag: v8.0.0~20488 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f3a356dbd499a7a1ddca0a2925d199911c0d7989;p=dealii.git scripts for assembling the restriction matrices for DG elements. Restriction for DG is not an interpolation but a projection git-svn-id: https://svn.dealii.org/trunk@2918 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/scripts/1d/restriction_dg b/deal.II/deal.II/source/fe/scripts/1d/restriction_dg new file mode 100644 index 0000000000..6bdcef8b25 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/1d/restriction_dg @@ -0,0 +1,100 @@ + 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: + + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/1d/restriction_dg1 b/deal.II/deal.II/source/fe/scripts/1d/restriction_dg1 new file mode 100644 index 0000000000..5e143fea43 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/1d/restriction_dg1 @@ -0,0 +1,40 @@ +# --------------------------------- 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); + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/1d/restriction_dg2 b/deal.II/deal.II/source/fe/scripts/1d/restriction_dg2 new file mode 100644 index 0000000000..75639c3018 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/1d/restriction_dg2 @@ -0,0 +1,41 @@ +# --------------------------------- 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); + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/1d/restriction_dg3 b/deal.II/deal.II/source/fe/scripts/1d/restriction_dg3 new file mode 100644 index 0000000000..ba3dd1745e --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/1d/restriction_dg3 @@ -0,0 +1,42 @@ +# --------------------------------- 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); + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/1d/restriction_dg4 b/deal.II/deal.II/source/fe/scripts/1d/restriction_dg4 new file mode 100644 index 0000000000..b3b2c34817 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/1d/restriction_dg4 @@ -0,0 +1,43 @@ +# --------------------------------- 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); + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/2d/restriction_dg b/deal.II/deal.II/source/fe/scripts/2d/restriction_dg new file mode 100644 index 0000000000..a789e31004 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/2d/restriction_dg @@ -0,0 +1,121 @@ + 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: + + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/2d/restriction_dg1 b/deal.II/deal.II/source/fe/scripts/2d/restriction_dg1 new file mode 100644 index 0000000000..bf8c1da2ca --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/2d/restriction_dg1 @@ -0,0 +1,44 @@ +# --------------------------------- 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); + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/2d/restriction_dg2 b/deal.II/deal.II/source/fe/scripts/2d/restriction_dg2 new file mode 100644 index 0000000000..c093e67d9b --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/2d/restriction_dg2 @@ -0,0 +1,40 @@ +# --------------------------------- 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); + diff --git a/deal.II/deal.II/source/fe/scripts/2d/restriction_dg3 b/deal.II/deal.II/source/fe/scripts/2d/restriction_dg3 new file mode 100644 index 0000000000..f4d60ec4ad --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/2d/restriction_dg3 @@ -0,0 +1,58 @@ +# --------------------------------- 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); + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/2d/restriction_dg4 b/deal.II/deal.II/source/fe/scripts/2d/restriction_dg4 new file mode 100644 index 0000000000..3f4480053e --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/2d/restriction_dg4 @@ -0,0 +1,68 @@ +# --------------------------------- 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); + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/3d/restriction_dg b/deal.II/deal.II/source/fe/scripts/3d/restriction_dg new file mode 100644 index 0000000000..9bcc76ec55 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/3d/restriction_dg @@ -0,0 +1,146 @@ + 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: + + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/3d/restriction_dg1 b/deal.II/deal.II/source/fe/scripts/3d/restriction_dg1 new file mode 100644 index 0000000000..05dae43b0c --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/3d/restriction_dg1 @@ -0,0 +1,50 @@ +# --------------------------------- 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); + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/3d/restriction_dg2 b/deal.II/deal.II/source/fe/scripts/3d/restriction_dg2 new file mode 100644 index 0000000000..14dcd93dcb --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/3d/restriction_dg2 @@ -0,0 +1,59 @@ +# --------------------------------- 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); + + + + + + + + + + + diff --git a/deal.II/deal.II/source/fe/scripts/3d/restriction_dg3 b/deal.II/deal.II/source/fe/scripts/3d/restriction_dg3 new file mode 100644 index 0000000000..a315aa5422 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/3d/restriction_dg3 @@ -0,0 +1,62 @@ +# --------------------------------- 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); + + + + + + + + + + +