From: hartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 22 May 2000 12:38:30 +0000 (+0000)
Subject: scripts for assembling the restriction matrices for DG elements. Restriction for... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e5fa56ce76f945b0e8933a54a51b2180b0be4df9;p=dealii-svn.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);
+
+
+
+
+
+
+
+
+
+
+