]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
First versions of these files.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 12 Feb 1999 09:24:26 +0000 (09:24 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 12 Feb 1999 09:24:26 +0000 (09:24 +0000)
git-svn-id: https://svn.dealii.org/trunk@785 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/scripts/3d/computations_on_faces [new file with mode: 0644]
deal.II/deal.II/source/fe/scripts/3d/lagrange [new file with mode: 0644]
deal.II/deal.II/source/fe/scripts/3d/lagrange-linear [new file with mode: 0644]
deal.II/deal.II/source/fe/scripts/3d/lagrange-quadratic [new file with mode: 0644]
deal.II/deal.II/source/fe/scripts/3d/lagrange-tools [new file with mode: 0644]
deal.II/deal.II/source/fe/scripts/3d/postprocess [new file with mode: 0644]

diff --git a/deal.II/deal.II/source/fe/scripts/3d/computations_on_faces b/deal.II/deal.II/source/fe/scripts/3d/computations_on_faces
new file mode 100644 (file)
index 0000000..019eece
--- /dev/null
@@ -0,0 +1,31 @@
+  # 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);
+  z := 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):
+  z_real := sum(z[s]*tphi[s], s=0..3):
+
+  image := vector([x_real, y_real, z_real]):
+
+  outward_vector := linalg[crossprod] (map(diff, image, xi), 
+                                       map(diff, image,eta)):
+  detJ := linalg[norm] (outward_vector, 2):
+  normal_vector := map (p->p/detJ, outward_vector):
+
+
+  readlib(C):
+  C(detJ, optimized):
+
+  # apply the following perl scripts:
+  # perl -pi -e 's/x\[(\d)\]/vertices[$1](0)/g; s/y\[(\d)\]/vertices[$1](1)/g; s/z\[(\d)\]/vertices[$1](2)/g;'
+  # perl -pi -e 's/^\s*t/const double t/g;'
\ No newline at end of file
diff --git a/deal.II/deal.II/source/fe/scripts/3d/lagrange b/deal.II/deal.II/source/fe/scripts/3d/lagrange
new file mode 100644 (file)
index 0000000..7dc2549
--- /dev/null
@@ -0,0 +1,248 @@
+# Maple script to compute much of the data needed to implement the
+# family of Lagrange elements in 3d. Expects that the fields denoting
+# position and number of support points, etc are already set. Note that
+# we assume a bilinear mapping from the unit to the real cell.
+#
+# $Id$
+# Author: Wolfgang Bangerth, 1998
+
+  phi_polynom := array(0..n_functions-1):
+  grad_phi_polynom := array(0..n_functions-1,0..2):
+  grad_grad_phi_polynom := array(0..n_functions-1,0..2,0..2):
+  local_mass_matrix := array(0..n_functions-1, 0..n_functions-1):
+  prolongation := array(0..7,0..n_functions-1, 0..n_functions-1):
+  interface_constraints := array(0..n_constraints-1,
+                                 0..n_face_functions-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],
+                  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):
+    grad_phi_polynom[i,0] := diff(phi_polynom[i], xi):
+    grad_phi_polynom[i,1] := diff(phi_polynom[i], eta):
+    grad_phi_polynom[i,2] := diff(phi_polynom[i], zeta):
+
+    grad_grad_phi_polynom[i,0,0] := diff(phi_polynom[i], xi, xi):
+    grad_grad_phi_polynom[i,0,1] := diff(phi_polynom[i], xi, eta):
+    grad_grad_phi_polynom[i,0,2] := diff(phi_polynom[i], xi, zeta):
+    grad_grad_phi_polynom[i,1,0] := diff(phi_polynom[i], eta,xi):
+    grad_grad_phi_polynom[i,1,1] := diff(phi_polynom[i], eta,eta):
+    grad_grad_phi_polynom[i,1,2] := diff(phi_polynom[i], eta,zeta):
+    grad_grad_phi_polynom[i,2,0] := diff(phi_polynom[i], zeta,xi):
+    grad_grad_phi_polynom[i,2,1] := diff(phi_polynom[i], zeta,eta):
+    grad_grad_phi_polynom[i,2,2] := diff(phi_polynom[i], zeta,zeta):
+  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
+  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:  
+
+  print ("Computing prolongation matrices"):
+  for i from 0 to 7 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], points[i][j,3]):
+      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, z)
+                    if ((x>1/2) or (y>1/2) or (z>1/2)) then
+                     0:
+                   else
+                     phi(i,2*x,2*y,2*z):
+                   fi:
+                 end: 
+  child_phi[1] := proc(i, x, y, z)
+                    if ((x<1/2) or (y>1/2) or (z>1/2)) then
+                     0:
+                   else
+                     phi(i,2*x-1,2*y, 2*z):
+                   fi:
+                 end: 
+  child_phi[2] := proc(i, x, y, z)
+                    if ((x<1/2) or (y>1/2) or (z<1/2)) then
+                     0:
+                   else
+                     phi(i,2*x-1,2*y, 2*z-1):
+                   fi:
+                 end: 
+  child_phi[3] := proc(i, x, y, z)
+                    if ((x>1/2) or (y>1/2) or (z<1/2)) then
+                     0:
+                   else
+                     phi(i,2*x,2*y,2*z-1):
+                   fi:
+                 end: 
+  child_phi[4] := proc(i, x, y, z)
+                    if ((x>1/2) or (y<1/2) or (z>1/2)) then
+                     0:
+                   else
+                     phi(i,2*x,2*y-1,2*z):
+                   fi:
+                 end: 
+  child_phi[5] := proc(i, x, y, z)
+                    if ((x<1/2) or (y<1/2) or (z>1/2)) then
+                     0:
+                   else
+                     phi(i,2*x-1,2*y-1,2*z):
+                   fi:
+                 end: 
+  child_phi[6] := proc(i, x, y, z)
+                    if ((x<1/2) or (y<1/2) or (z<1/2)) then
+                     0:
+                   else
+                     phi(i,2*x-1,2*y-1,2*z-1):
+                   fi:
+                 end: 
+  child_phi[7] := proc(i, x, y, z)
+                    if ((x>1/2) or (y<1/2) or (z<1/2)) then
+                     0:
+                   else
+                     phi(i,2*x,2*y-1,2*z-1):
+                   fi:
+                 end: 
+  restriction := array(0..7,0..n_functions-1, 0..n_functions-1):
+  for child from 0 to 7 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],
+                                                  support_points[j][3]):
+      od:
+    od:
+  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
+  #
+  # first compute for each function on the (large) face the polynom
+  # we get this by evaluating the respective global trial function
+  # with y=0
+  face_phi_polynom := array(0..n_face_functions-1):
+  for j from 0 to n_face_functions-1 do
+    face_phi_polynom[j] := proc(xi,eta)
+      subs(dummy=0, phi(constrained_face_function[j],xi,dummy,eta)):
+    end:
+  od:
+
+  for i from 0 to n_constraints-1 do
+    for j from 0 to n_face_functions-1 do
+      interface_constraints[i,j] 
+            := face_phi_polynom[j](constrained_face_support_points[i][0],
+                                   constrained_face_support_points[i][1]):
+    od:
+  od:
+
+                                     
+  # 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. same for z
+  #
+  print ("Computing real space support points"):
+  x := array(0..7);
+  y := array(0..7);
+  z := array(0..7):
+  tphi[0] := (1-xi)*(1-eta)*(1-zeta):
+  tphi[1] := xi*(1-eta)*(1-zeta):
+  tphi[2] := xi*(1-eta)*zeta:
+  tphi[3] := (1-xi)*(1-eta)*zeta:
+  tphi[4] := (1-xi)*eta*(1-zeta):
+  tphi[5] := xi*eta*(1-zeta):
+  tphi[6] := xi*eta*zeta:
+  tphi[7] := (1-xi)*eta*zeta:
+  x_real := sum(x[s]*tphi[s], s=0..7):
+  y_real := sum(y[s]*tphi[s], s=0..7):
+  z_real := sum(z[s]*tphi[s], s=0..7):
+
+  real_space_points := array(0..n_functions-1,0..2):
+  for i from 0 to n_functions-1 do
+    real_space_points[i,0] :=
+                   subs(xi=support_points[i][1],
+                        eta=support_points[i][2],
+                        zeta=support_points[i][3],
+                        x_real):
+    real_space_points[i,1] :=
+                   subs(xi=support_points[i][1],
+                        eta=support_points[i][2],
+                        zeta=support_points[i][3],
+                        y_real):
+    real_space_points[i,2] :=
+                   subs(xi=support_points[i][1],
+                        eta=support_points[i][2],
+                        zeta=support_points[i][3],
+                        z_real):
+   od:
\ No newline at end of file
diff --git a/deal.II/deal.II/source/fe/scripts/3d/lagrange-linear b/deal.II/deal.II/source/fe/scripts/3d/lagrange-linear
new file mode 100644 (file)
index 0000000..f33f025
--- /dev/null
@@ -0,0 +1,125 @@
+#  --------------------------------- For 3d ---------------------------------
+#  -- Use the following maple script to generate the basis functions,
+#  -- gradients and prolongation matrices as well as the mass matrix.
+#  -- Make sure that the files do not exists beforehand, since output
+#  -- is appended instead of overwriting previous contents.
+#  --
+#  -- You should only have to change the very first lines for polynomials
+#  -- of higher order.
+#
+#  --------------------------------------------------------------------------
+#
+# $Id$
+# Author: Wolfgang Bangerth, 1999
+
+  n_functions      := 8:
+  n_face_functions := 4:
+  n_constraints    := 5:
+
+  trial_function := ((a1 + a2*xi) +
+                     (b1 + b2*xi)*eta) +
+                    ((d1 + d2*xi) + 
+                     (e1 + e2*xi)*eta)*zeta:
+  face_trial_function := subs(zeta=0, trial_function):
+  # note: support_points[i] is a vector which is indexed from
+  # one and not from zero!
+  support_points := array(0..n_functions-1):
+  support_points[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]):
+
+  face_support_points := array(0..n_face_functions-1):
+  face_support_points[0] := [0,0]:
+  face_support_points[1] := [1,0]:
+  face_support_points[2] := [1,1]:
+  face_support_points[3] := [0,1]:
+
+  # list of functions which are at face 0, used to compute
+  # the constraints on a face
+  constrained_face_function := array (0..n_face_functions-1):
+  # the list of points at which we want the functions at
+  # faces to be evaluated
+  constrained_face_support_points := array(0..n_constraints-1):
+  constrained_face_function[0] := 0:
+  constrained_face_function[1] := 1:
+  constrained_face_function[2] := 2:
+  constrained_face_function[3] := 3:
+  constrained_face_support_points[0] := array(0..1, [1/2,1/2]):
+  constrained_face_support_points[1] := array(0..1, [1/2,0]):
+  constrained_face_support_points[2] := array(0..1, [1,1/2]):
+  constrained_face_support_points[3] := array(0..1, [1/2,1]):
+  constrained_face_support_points[4] := array(0..1, [0,1/2]):
+
+
+  # do the real work
+  read "lagrange":
+
+  # ... originally taken from another comment, so this does not
+  # fit in here too well...
+  #
+  # Since we're already at it and need it anyway, we also compute the
+  # Jacobian matrix of the transform and its derivatives. For the
+  # question of whether to take the given form or its transpose, refer
+  # to the documentation of the FEValues class and the source code
+  # documentation of FELinearMapping::fill_fe_values. Also note, that
+  # the computed inverse is multiplied to the unit cell gradients
+  # *from the right*. 
+  print ("Computing Jacobian matrices"):
+  Jacobian := linalg[matrix](3,3, [[diff(x_real,xi), diff(x_real,eta), diff(x_real,zeta)],
+                                   [diff(y_real,xi), diff(y_real,eta), diff(y_real,zeta)],
+                                   [diff(z_real,xi), diff(z_real,eta), diff(z_real,zeta)]]):
+  inverseJacobian := linalg[inverse](Jacobian):
+  detJ := linalg[det](Jacobian):
+
+  grad_inverseJacobian := array(1..3, 1..3, 1..3):
+  for i from 1 to 3 do
+    for j from 1 to 3 do
+      for k from 1 to 3 do
+        if (i=1) then
+          grad_inverseJacobian[i,j,k] := diff (inverseJacobian[j,k], xi):
+        else 
+          if (i=2) then
+            grad_inverseJacobian[i,j,k] := diff (inverseJacobian[j,k], eta):
+          else
+            grad_inverseJacobian[i,j,k] := diff (inverseJacobian[j,k], zeta):
+         fi:
+        fi:
+      od:
+    od:
+  od:
+
+
+  print ("computing support points in real space"):
+  real_points := array(0..n_functions-1, 0..2);
+  for i from 0 to n_functions-1 do
+    real_points[i,0] := subs(xi=support_points[i][1],
+                            eta=support_points[i][2],  
+                             zeta=support_points[i][3], x_real);
+    real_points[i,1] := subs(xi=support_points[i][1],
+                             eta=support_points[i][2], 
+                             zeta=support_points[i][3], y_real);
+    real_points[i,2] := subs(xi=support_points[i][1],
+                             eta=support_points[i][2], 
+                             zeta=support_points[i][3], z_real);
+  od:
+
+
+  
+  # write data to files
+  print ("writing data to files"):
+  readlib(C):
+  C(phi_polynom,            filename=linear3d_shape_value):
+  C(grad_phi_polynom,       filename=linear3d_shape_grad):
+  C(grad_grad_phi_polynom,  filename=linear3d_shape_grad_grad):
+  C(prolongation,           filename=linear3d_prolongation):
+  C(restriction,            filename=linear3d_restriction):
+  C(local_mass_matrix,      filename=linear3d_massmatrix):
+  C(interface_constraints,  filename=linear3d_constraints):
+  C(real_points, optimized, filename=linear3d_real_points):
+  C(inverseJacobian,        filename=linear3d_inverse_jacobian):
+  C(grad_inverseJacobian,   filename=linear3d_inverse_jacobian_grad):
diff --git a/deal.II/deal.II/source/fe/scripts/3d/lagrange-quadratic b/deal.II/deal.II/source/fe/scripts/3d/lagrange-quadratic
new file mode 100644 (file)
index 0000000..69dfc0a
--- /dev/null
@@ -0,0 +1,113 @@
+#  --------------------------------- For 3d ---------------------------------
+#  -- Use the following maple script to generate the basis functions,
+#  -- gradients and prolongation matrices as well as the mass matrix.
+#  -- Make sure that the files do not exists beforehand, since output
+#  -- is appended instead of overwriting previous contents.
+#  --
+#  -- You should only have to change the very first lines for polynomials
+#  -- of higher order.
+#
+#  --------------------------------------------------------------------------
+#
+# $Id$
+# Author: Wolfgang Bangerth, 1999
+
+  read "lagrange-tools":
+
+  n_functions      := 27:
+  n_face_functions := 9:
+  n_constraints    := 21:
+
+  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:
+  face_trial_function := subs(zeta=0, trial_function):
+  # note: support_points[i] is a vector which is indexed from
+  # one and not from zero!
+  support_points := array(0..n_functions-1):
+
+  support_points_fill_vertices (0, support_points):
+  support_points_fill_lines (8, 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
+
+  face_support_points := array(0..n_face_functions-1):
+  face_support_points[0] := [0,0]:
+  face_support_points[1] := [1,0]:
+  face_support_points[2] := [1,1]:
+  face_support_points[3] := [0,1]:
+  face_support_points[4] := [1/2,0]:
+  face_support_points[5] := [1,1/2]:
+  face_support_points[6] := [1/2,1]:
+  face_support_points[7] := [0,1/2]:
+  face_support_points[8] := [1/2,1/2]:
+
+  # list of functions which are at face 0, used to compute
+  # the constraints on a face
+  constrained_face_function := array (0..n_face_functions-1):
+  # the list of points at which we want the functions at
+  # faces to be evaluated
+  constrained_face_support_points := array(0..n_constraints-1):
+  constrained_face_function[0] := 0:
+  constrained_face_function[1] := 1:
+  constrained_face_function[2] := 2:
+  constrained_face_function[3] := 3:
+  constrained_face_function[4] := 8:
+  constrained_face_function[5] := 9:
+  constrained_face_function[6] := 10:
+  constrained_face_function[7] := 11:
+  constrained_face_function[8] := 20:
+  constrained_face_support_points[0] := array(0..1, [1/2,1/2]):  # center vertex
+  constrained_face_support_points[1] := array(0..1, [1/2,0]):    # centers of large lines
+  constrained_face_support_points[2] := array(0..1, [1,1/2]):
+  constrained_face_support_points[3] := array(0..1, [1/2,1]):
+  constrained_face_support_points[4] := array(0..1, [0,1/2]):
+  constrained_face_support_points[5] := array(0..1, [1/2,1/4]):  # lines from center to boundary
+  constrained_face_support_points[6] := array(0..1, [3/4,1/2]):
+  constrained_face_support_points[7] := array(0..1, [1/2,3/4]):
+  constrained_face_support_points[8] := array(0..1, [1/4,1/2]):
+  constrained_face_support_points[9] := array(0..1, [1/4,0]):    # children of bounding lines
+  constrained_face_support_points[10] := array(0..1, [3/4,0]):
+  constrained_face_support_points[11] := array(0..1, [1,1/4]):
+  constrained_face_support_points[12] := array(0..1, [1,3/4]):
+  constrained_face_support_points[13] := array(0..1, [1/4,1]):
+  constrained_face_support_points[14] := array(0..1, [3/4,1]):
+  constrained_face_support_points[15] := array(0..1, [0,1/4]):
+  constrained_face_support_points[16] := array(0..1, [0,3/4]):
+  constrained_face_support_points[17] := array(0..1, [1/4,1/4]):    # child quads
+  constrained_face_support_points[18] := array(0..1, [3/4,1/4]):
+  constrained_face_support_points[19] := array(0..1, [3/4,3/4]):
+  constrained_face_support_points[20] := array(0..1, [1/4,3/4]):
+
+
+  # do the real work
+  read "lagrange":
+
+
+  
+  # write data to files
+  print ("writing data to files"):
+  readlib(C):
+  C(phi_polynom,            filename=quadratic3d_shape_value):
+  C(grad_phi_polynom,       filename=quadratic3d_shape_grad):
+  C(grad_grad_phi_polynom,  filename=quadratic3d_shape_grad_grad):
+  C(prolongation,           filename=quadratic3d_prolongation):
+  C(restriction,            filename=quadratic3d_restriction):
+  C(interface_constraints,  filename=quadratic3d_constraints):
+  C(real_space_points,      optimized, filename=quadratic3d_real_points):
+
+  writeto (quadratic3d_unit_support_points):
+  print (support_points):
+
+  
\ No newline at end of file
diff --git a/deal.II/deal.II/source/fe/scripts/3d/lagrange-tools b/deal.II/deal.II/source/fe/scripts/3d/lagrange-tools
new file mode 100644 (file)
index 0000000..219b068
--- /dev/null
@@ -0,0 +1,106 @@
+support_points_fill_vertices := proc (starting_index, support_points)
+  support_points[starting_index+0] := array(1..3, [0,0,0]):
+  support_points[starting_index+1] := array(1..3, [1,0,0]):
+  support_points[starting_index+2] := array(1..3, [1,0,1]):
+  support_points[starting_index+3] := array(1..3, [0,0,1]):
+  support_points[starting_index+4] := array(1..3, [0,1,0]):
+  support_points[starting_index+5] := array(1..3, [1,1,0]):
+  support_points[starting_index+6] := array(1..3, [1,1,1]):
+  support_points[starting_index+7] := array(1..3, [0,1,1]):
+end:
+
+
+
+support_points_fill_lines := proc (starting_index, dofs_per_line, support_points)
+
+  local next_index, increment, i:
+
+  next_index := starting_index:
+  increment := 1/(dofs_per_line+1):
+
+  # line 0
+  for i from 1 to dofs_per_line do
+    print (next_index+i-1):
+    support_points[next_index+i-1] 
+      := array (1..3, [i*increment, 0, 0]):
+    next_index := next_index+1
+  od:
+
+  # line 1
+  for i from 1 to dofs_per_line do
+    support_points[next_index+i-1] 
+      := array (1..3, [1, 0, i*increment]):
+    next_index := next_index+1
+  od:
+
+  # line 2
+  for i from 1 to dofs_per_line do
+    support_points[next_index+i-1] 
+      := array (1..3, [i*increment, 0, 1]):
+    next_index := next_index+1
+  od:
+
+  # line 3
+  for i from 1 to dofs_per_line do
+    support_points[next_index+i-1] 
+      := array (1..3, [0, 0, i*increment]):
+    next_index := next_index+1
+  od:
+
+  # line 4
+  for i from 1 to dofs_per_line do
+    support_points[next_index+i-1] 
+      := array (1..3, [i*increment, 1, 0]):
+    next_index := next_index+1
+  od:
+
+  # line 5
+  for i from 1 to dofs_per_line do
+    support_points[next_index+i-1] 
+      := array (1..3, [1, 1, i*increment]):
+    next_index := next_index+1
+  od:
+
+  # line 6
+  for i from 1 to dofs_per_line do
+    support_points[next_index+i-1] 
+      := array (1..3, [i*increment, 1, 1]):
+    next_index := next_index+1
+  od:
+
+  # line 7
+  for i from 1 to dofs_per_line do
+    support_points[next_index+i-1] 
+      := array (1..3, [0, 1, i*increment]):
+    next_index := next_index+1
+  od:
+
+
+  # line 8
+  for i from 1 to dofs_per_line do
+    support_points[next_index+i-1] 
+      := array (1..3, [0, i*increment,0]):
+    next_index := next_index+1
+  od:
+
+  # line 9
+  for i from 1 to dofs_per_line do
+    support_points[next_index+i-1] 
+      := array (1..3, [1, i*increment, 0]):
+    next_index := next_index+1
+  od:
+
+  # line 10
+  for i from 1 to dofs_per_line do
+    support_points[next_index+i-1] 
+      := array (1..3, [1, i*increment, 1]):
+    next_index := next_index+1
+  od:
+
+  # line 11
+  for i from 1 to dofs_per_line do
+    support_points[next_index+i-1] 
+      := array (1..3, [0, i*increment, 1]):
+    next_index := next_index+1
+  od:
+end:
\ No newline at end of file
diff --git a/deal.II/deal.II/source/fe/scripts/3d/postprocess b/deal.II/deal.II/source/fe/scripts/3d/postprocess
new file mode 100644 (file)
index 0000000..d783b63
--- /dev/null
@@ -0,0 +1,56 @@
+#  Use the following perl scripts to convert the output into a
+#  suitable format.
+#
+#  $Id$
+#  Wolfgang Bangerth, 1998
+  
+perl -pi -e 's/phi_polynom\[(\d+)\] =/case $1: return/g;' *3d_shape_value
+
+perl -pi -e 's/([^;])\n/$1/g;' *3d_shape_grad
+perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[0\] = (.*);/case $1: return Point<3>($2,/g;' *3d_shape_grad
+perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[[01]\] = (.*);/$2,/g;' *3d_shape_grad
+perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[2\] = (.*);/$2);/g;' *3d_shape_grad
+
+
+#    concatenate all lines for each entry
+perl -pi -e 's/([^;])\n/$1/g;' *3d_shape_grad_grad
+#    rename the variable
+perl -pi -e 's/\s*grad_grad_phi_polynom/return_value/g;' *3d_shape_grad_grad
+#    insert 'case' and 'break' statements
+perl -pi -e 's/(return_value\[(\d)\]\[0\]\[0\] = .*;)/break;\ncase $2:\n$1/g;' *3d_shape_grad_grad
+#    eliminate first index, since that one is caught by the 'case' statement
+perl -pi -e 's/return_value\[\d+\]/return_value/g;' *3d_shape_grad_grad
+#    delete lines where only a zero is set, since this already is done in the constructor
+perl -pi -e 's/.*= 0.0;\n//g;' *3d_shape_grad_grad
+
+
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' *3d_prolongation
+perl -pi -e 's/.*= 0.0;\n//g;' *3d_prolongation
+
+
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' *3d_restriction
+perl -pi -e 's/.*= 0.0;\n//g;' *3d_restriction
+
+
+perl -pi -e 's/\[(\d+)\]\[(\d+)\]/($1,$2)/g;' *3d_constraints
+perl -pi -e 's/.*= 0.0;\n//g;' *3d_constraints
+
+perl -pi -e 's/^\s*t/const double t/g;' *3d_inverse_jacobian
+perl -pi -e 's/x\[(\d)\]/vertices[$1](0)/g;' *3d_inverse_jacobian
+perl -pi -e 's/y\[(\d)\]/vertices[$1](1)/g;' *3d_inverse_jacobian
+perl -pi -e 's/inverseJacobian/jacobians[point]/g;' *3d_inverse_jacobian
+perl -pi -e 's/\[(\d)\]\[(\d)\] =/($1,$2) =/g;' *3d_inverse_jacobian
+
+perl -pi -e 's/^\s*t/const double t/g;' *3d_inverse_jacobian_grad
+perl -pi -e 's/x\[(\d)\]/vertices[$1](0)/g;' *3d_inverse_jacobian_grad
+perl -pi -e 's/y\[(\d)\]/vertices[$1](1)/g;' *3d_inverse_jacobian_grad
+perl -pi -e 's/inverseJacobian/jacobians_grad[point]/g;' *3d_inverse_jacobian_grad
+
+
+perl -pi -e 's/^array.*\n//g; s/^\s*\]\)//g; s/^\n//g;' *3d_unit_support_points
+perl -pi -e 's/\s+\((\d+)\)/  unit_points[$1]/g;' *3d_unit_support_points
+perl -pi -e 's/= \[/= Point<3>(/g; s/\]\s*\n/);\n/g;' *3d_unit_support_points
+
+
+perl -pi -e 's/real_space_points\[(\d+)\]\[(\d+)\]/support_points[$1]($2)/g;' *3d_real_points
+perl -pi -e 's/x\[(\d+)\]/vertices[$1](0)/g; s/y\[(\d+)\]/vertices[$1](1)/g; s/z\[(\d+)\]/vertices[$1](2)/g;' *3d_real_points
\ No newline at end of file

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.