From 0d23044fcc339f00c609c445d3659e99d9446051 Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 4 Nov 1998 10:38:43 +0000 Subject: [PATCH] Separate the maple scripts to generate the finite elements from the .cc files. The reason is that much of the code is the same for all elements and when the code grew larger, maintenance of these parts grew more and more difficult. Now, we use the maple 'read' statement and single out the common stuff into the file 'lagrange'. git-svn-id: https://svn.dealii.org/trunk@629 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/scripts/2d/lagrange | 177 ++++++++++++++++++ .../source/fe/scripts/2d/lagrange-cubic | 68 +++++++ .../source/fe/scripts/2d/lagrange-quadratic | 66 +++++++ .../source/fe/scripts/2d/lagrange-quartic | 82 ++++++++ .../deal.II/source/fe/scripts/2d/postprocess | 14 ++ 5 files changed, 407 insertions(+) create mode 100644 deal.II/deal.II/source/fe/scripts/2d/lagrange create mode 100644 deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic create mode 100644 deal.II/deal.II/source/fe/scripts/2d/lagrange-quadratic create mode 100644 deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic create mode 100644 deal.II/deal.II/source/fe/scripts/2d/postprocess 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 index 0000000000..b66816fe0e --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/2d/lagrange @@ -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 index 0000000000..431b6781c7 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic @@ -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 index 0000000000..deb5404622 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/2d/lagrange-quadratic @@ -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 index 0000000000..d6ced17de0 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic @@ -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 index 0000000000..a19e779914 --- /dev/null +++ b/deal.II/deal.II/source/fe/scripts/2d/postprocess @@ -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 -- 2.39.5