From 2beeb9d4019e484cdadf5f169c342c0ee2b36dc4 Mon Sep 17 00:00:00 2001 From: wolf Date: Fri, 12 Feb 1999 09:24:26 +0000 Subject: [PATCH] First versions of these files. git-svn-id: https://svn.dealii.org/trunk@785 0785d39b-7218-0410-832d-ea1e28bc413d --- .../fe/scripts/3d/computations_on_faces | 31 +++ deal.II/deal.II/source/fe/scripts/3d/lagrange | 248 ++++++++++++++++++ .../source/fe/scripts/3d/lagrange-linear | 125 +++++++++ .../source/fe/scripts/3d/lagrange-quadratic | 113 ++++++++ .../source/fe/scripts/3d/lagrange-tools | 106 ++++++++ .../deal.II/source/fe/scripts/3d/postprocess | 56 ++++ 6 files changed, 679 insertions(+) create mode 100644 deal.II/deal.II/source/fe/scripts/3d/computations_on_faces create mode 100644 deal.II/deal.II/source/fe/scripts/3d/lagrange create mode 100644 deal.II/deal.II/source/fe/scripts/3d/lagrange-linear create mode 100644 deal.II/deal.II/source/fe/scripts/3d/lagrange-quadratic create mode 100644 deal.II/deal.II/source/fe/scripts/3d/lagrange-tools create mode 100644 deal.II/deal.II/source/fe/scripts/3d/postprocess 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 index 0000000000..019eece09a --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/3d/computations_on_faces @@ -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 index 0000000000..7dc2549e7b --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/3d/lagrange @@ -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 index 0000000000..f33f0257aa --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/3d/lagrange-linear @@ -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 index 0000000000..69dfc0adda --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/3d/lagrange-quadratic @@ -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 index 0000000000..219b068245 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/3d/lagrange-tools @@ -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 index 0000000000..d783b63f0e --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/3d/postprocess @@ -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 -- 2.39.5