]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Separate the maple scripts to generate the finite elements from the .cc files. The...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 4 Nov 1998 10:38:43 +0000 (10:38 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 4 Nov 1998 10:38:43 +0000 (10:38 +0000)
git-svn-id: https://svn.dealii.org/trunk@629 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/scripts/2d/lagrange [new file with mode: 0644]
deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic [new file with mode: 0644]
deal.II/deal.II/source/fe/scripts/2d/lagrange-quadratic [new file with mode: 0644]
deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic [new file with mode: 0644]
deal.II/deal.II/source/fe/scripts/2d/postprocess [new file with mode: 0644]

diff --git a/deal.II/deal.II/source/fe/scripts/2d/lagrange b/deal.II/deal.II/source/fe/scripts/2d/lagrange
new file mode 100644 (file)
index 0000000..b66816f
--- /dev/null
@@ -0,0 +1,177 @@
+# Maple script to compute much of the data needed to implement the
+# family of Lagrange elements in 2d. Expects that the fields denoting
+# position and number of support points, etc are already set.
+#
+# $Id$
+# Author: Wolfgang Bangerth, 1998
+
+  phi_polynom := array(0..n_functions-1):
+  grad_phi_polynom := array(0..n_functions-1,0..1):
+  local_mass_matrix := array(0..n_functions-1, 0..n_functions-1):
+  prolongation := array(0..3,0..n_functions-1, 0..n_functions-1):
+  interface_constraints := array(0..2*(n_face_functions-2)+1-1,
+                                 0..n_face_functions-1):
+  real_points := array(0..n_functions-1, 0..1);
+
+  print ("Computing basis functions"):
+  for i from 0 to n_functions-1 do
+    print (i):
+    values := array(1..n_functions):
+    for j from 1 to n_functions do
+      values[j] := 0:
+    od:  
+    values[i+1] := 1:
+
+    equation_system := {}:
+    for j from 0 to n_functions-1 do
+      poly := subs(xi=support_points[j][1],
+                   eta=support_points[j][2],
+                  trial_function):
+      if (i=j) then
+        equation_system := equation_system union {poly = 1}:
+      else     
+        equation_system := equation_system union {poly = 0}:
+      fi:      
+    od:
+    
+    phi_polynom[i] := subs(solve(equation_system), trial_function):
+    grad_phi_polynom[i,0] := diff(phi_polynom[i], xi):
+    grad_phi_polynom[i,1] := diff(phi_polynom[i], eta):
+  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
+  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:  
+
+  print ("Computing prolongation matrices"):
+  for i from 0 to 3 do
+    print ("child", i):
+    for j from 0 to n_functions-1 do
+      for k from 0 to n_functions-1 do
+        prolongation[i,j,k] := phi(k, points[i][j,1], points[i][j,2]):
+      od:
+    od:
+  od:
+
+  print ("Computing restriction matrices"):
+  # to get the restriction (interpolation) matrices, evaluate
+  # the basis functions on the child cells at the global
+  # interpolation points
+  child_phi[0] := proc(i, x, y)
+                    if ((x>1/2) or (y>1/2)) then
+                     0:
+                   else
+                     phi(i,2*x,2*y):
+                   fi:
+                 end: 
+  child_phi[1] := proc(i, x, y)
+                    if ((x<1/2) or (y>1/2)) then
+                     0:
+                   else
+                     phi(i,2*x-1,2*y):
+                   fi:
+                 end: 
+  child_phi[2] := proc(i, x, y)
+                    if ((x<1/2) or (y<1/2)) then
+                     0:
+                   else
+                     phi(i,2*x-1,2*y-1):
+                   fi:
+                 end: 
+  child_phi[3] := proc(i, x, y)
+                    if ((x>1/2) or (y<1/2)) then
+                     0:
+                   else
+                     phi(i,2*x,2*y-1):
+                   fi:
+                 end: 
+  restriction := array(0..3,0..n_functions-1, 0..n_functions-1):
+  for child from 0 to 3 do
+    for j from 0 to n_functions-1 do
+      for k from 0 to n_functions-1 do
+        restriction[child,j,k] := child_phi[child](k,
+                                                  support_points[j][1],
+                                                  support_points[j][2]):
+      od:
+    od:
+  od:
+
+
+  print ("Computing local mass matrix"):
+  # tphi are the basis functions of the linear element. These functions
+  # are used for the computation of the subparametric transformation from
+  # unit cell to real cell.
+  # x and y are arrays holding the x- and y-values of the four vertices
+  # of this cell in real space. 
+  x := array(0..3);
+  y := array(0..3);
+  tphi[0] := (1-xi)*(1-eta):
+  tphi[1] := xi*(1-eta):
+  tphi[2] := xi*eta:
+  tphi[3] := (1-xi)*eta:
+  x_real := sum(x[s]*tphi[s], s=0..3):
+  y_real := sum(y[s]*tphi[s], s=0..3):
+  detJ := diff(x_real,xi)*diff(y_real,eta) - diff(x_real,eta)*diff(y_real,xi):
+  for i from 0 to n_functions-1 do
+    print ("line", i):
+    for j from 0 to n_functions-1 do
+      local_mass_matrix[i,j] := int(int(phi_polynom[i] * phi_polynom[j] * detJ,
+                                        xi=0..1), eta=0..1):
+    od:
+  od:
+
+  print ("computing support points in real space"):
+  for i from 0 to n_functions-1 do
+    real_points[i,0] := subs(xi=support_points[i][1],
+                             eta=support_points[i][2], x_real);
+    real_points[i,1] := subs(xi=support_points[i][1],
+                             eta=support_points[i][2], y_real);
+  od:
+
+  print ("computing interface constraint matrices"):
+  # compute the interface constraint matrices. these are the values of the
+  # basis functions on the coarse cell's face at the points of the child
+  # cell's basis functions on the child faces
+  face_phi_polynom := array(0..n_face_functions-1):
+  for i from 0 to n_face_functions-1 do
+    # note that the interp function wants vectors indexed from
+    #   one and not from zero. 
+    values := array(1..n_face_functions):
+    for j from 1 to n_face_functions do
+      values[j] := 0:
+    od:  
+    values[i+1] := 1:
+
+    shifted_face_support_points := array (1..n_face_functions):
+    for j from 1 to n_face_functions do
+      shifted_face_support_points[j] := face_support_points[j-1]:
+    od:
+    
+    face_phi_polynom[i] := interp (shifted_face_support_points, values, xi):
+  od:
+
+  for i from 0 to 2*(n_face_functions-2)+1-1 do
+    for j from 0 to n_face_functions-1 do
+      interface_constraints[i,j] := subs(xi=constrained_face_support_points[i],
+                                     face_phi_polynom[j]); 
+    od:
+  od:
diff --git a/deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic b/deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic
new file mode 100644 (file)
index 0000000..431b678
--- /dev/null
@@ -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: Wolfgang Bangerth, 1998
+
+  n_functions      := 16:
+  n_face_functions := 4:
+
+  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:
+  face_trial_function := a + b*xi + c*xi*xi + d*xi*xi*xi:
+  # note: support_points[i] is a vector which is indexed from
+  # one and not from zero!
+  support_points := array(0..n_functions-1):
+  support_points[0] := [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]:
+
+  face_support_points := array(0..n_face_functions-1):
+  face_support_points[0] := 0:
+  face_support_points[1] := 1:
+  face_support_points[2] := 1/3:
+  face_support_points[3] := 2/3:
+  constrained_face_support_points := array(0..2*(n_face_functions-2)+1-1):
+  constrained_face_support_points[0] := 1/2:
+  constrained_face_support_points[1] := 1/6:
+  constrained_face_support_points[2] := 2/6:
+  constrained_face_support_points[3] := 4/6:
+  constrained_face_support_points[4] := 5/6:
+
+  # do the real work
+  read "lagrange"
+
+  
+  # write data to files
+  print ("writing data to files"):
+  readlib(C):
+  C(phi_polynom, filename=cubic2d.shape_value):
+  C(grad_phi_polynom, filename=cubic2d.shape_grad):
+  C(prolongation, filename=cubic2d.prolongation):
+  C(restriction, filename=cubic2d.restriction):
+  C(local_mass_matrix, optimized, filename=cubic2d.massmatrix):
+  C(interface_constraints, filename=cubic2d.constraints):
+  C(real_points, optimized, filename=cubic2d.real_points):
+
diff --git a/deal.II/deal.II/source/fe/scripts/2d/lagrange-quadratic b/deal.II/deal.II/source/fe/scripts/2d/lagrange-quadratic
new file mode 100644 (file)
index 0000000..deb5404
--- /dev/null
@@ -0,0 +1,66 @@
+#  --------------------------------- 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.
+#  --
+#  -- Please note:
+#  -- Apart from the restriction matrices, I did not initially use it; it is
+#  -- an adaption of the script for cubic and quartic elements. For
+#  -- some of the data, however, a smaller script is given in the
+#  -- FEQuadratic<2> constructor.
+#
+#  --------------------------------------------------------------------------
+#
+# $Id$
+# Author: Wolfgang Bangerth, 1998
+
+  n_functions      := 9:
+  n_face_functions := 3:
+
+  trial_function := (a1 + a2*xi + a3*xi*xi) +
+                     (b1 + b2*xi + b3*xi*xi)*eta +
+                    (c1 + c2*xi + c3*xi*xi)*eta*eta:
+  face_trial_function := a + b*xi + c*xi*xi:
+  # note: support_points[i] is a vector which is indexed from
+  # one and not from zero!
+  support_points := array(0..n_functions-1):
+  support_points[0] := [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]:
+
+  face_support_points := array(0..n_face_functions-1):
+  face_support_points[0] := 0:
+  face_support_points[1] := 1:
+  face_support_points[2] := 1/2:
+
+  constrained_face_support_points := array(0..2*(n_face_functions-2)+1-1):
+  constrained_face_support_points[0] := 1/2:
+  constrained_face_support_points[1] := 1/4:
+  constrained_face_support_points[2] := 3/4:
+
+
+  # do the real work
+  read "lagrange"
+
+  
+  # write data to files
+  print ("writing data to files"):
+  readlib(C):
+  C(phi_polynom, filename=quadratic2d.shape_value):
+  C(grad_phi_polynom, filename=quadratic2d.shape_grad):
+  C(prolongation, filename=quadratic2d.prolongation):
+  C(restriction, filename=quadratic2d.restriction):
+  C(local_mass_matrix, optimized, filename=quadratic2d.massmatrix):
+  C(interface_constraints, filename=quadratic2d.constraints):
+  C(real_points, optimized, filename=quadratic2d.real_points):
+
diff --git a/deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic b/deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic
new file mode 100644 (file)
index 0000000..d6ced17
--- /dev/null
@@ -0,0 +1,82 @@
+#  --------------------------------- 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: Wolfgang Bangerth, 1998
+
+  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:
+  face_trial_function := a + b*xi + c*xi*xi + d*xi**3 + e*xi**4:
+  # note: support_points[i] is a vector which is indexed from
+  # one and not from zero!
+  support_points := array(0..n_functions-1):
+  support_points[0] := [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]:
+
+  face_support_points := array(0..n_face_functions-1):
+  face_support_points[0] := 0:
+  face_support_points[1] := 1:
+  face_support_points[2] := 1/4:
+  face_support_points[3] := 2/4:
+  face_support_points[4] := 3/4:
+  constrained_face_support_points := array(0..2*(n_face_functions-2)+1-1):
+  constrained_face_support_points[0] := 1/2:
+  constrained_face_support_points[1] := 1/8:
+  constrained_face_support_points[2] := 2/8:
+  constrained_face_support_points[3] := 3/8:
+  constrained_face_support_points[4] := 5/8:
+  constrained_face_support_points[5] := 6/8:
+  constrained_face_support_points[6] := 7/8:
+
+
+  # do the real work
+  read "lagrange"
+
+  
+  # write data to files
+  print ("writing data to files"):
+  readlib(C):
+  C(phi_polynom, filename=quartic2d.shape_value):
+  C(grad_phi_polynom, filename=quartic2d.shape_grad):
+  C(prolongation, filename=quartic2d.prolongation):
+  C(restriction, filename=quartic2d.restriction):
+  C(local_mass_matrix, optimized, filename=quartic2d.massmatrix):
+  C(interface_constraints, filename=quartic2d.constraints):
+  C(real_points, optimized, filename=quartic2d.real_points):
+
diff --git a/deal.II/deal.II/source/fe/scripts/2d/postprocess b/deal.II/deal.II/source/fe/scripts/2d/postprocess
new file mode 100644 (file)
index 0000000..a19e779
--- /dev/null
@@ -0,0 +1,14 @@
+#  Use the following perl scripts to convert the output into a
+#  suitable format.
+  
+perl -pi -e 's/phi_polynom\[(\d+)\] =/case $1: return/g;' *2d.shape_value
+perl -pi -e 's/([^;])\n/$1/g;' *2d.shape_grad
+perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[0\] = (.*);/case $1: return Point<2>($2,/g;' *2d.shape_grad
+perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[1\] = (.*);/$2);/g;' *2d.shape_grad
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]/($1,$2)/g;' *2d.massmatrix
+perl -pi -e 's/(t\d+) =/const double $1 =/g;' *2d.massmatrix
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' *2d.prolongation
+perl -pi -e 's/.*= 0.0;\n//g;' *2d.prolongation
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' *2d.restriction
+perl -pi -e 's/.*= 0.0;\n//g;' *2d.restriction
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]/($1,$2)/g;' *2d.constraints

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.