From d791b415c87975d2a643268f64a7433cfca94d8d Mon Sep 17 00:00:00 2001 From: Ralf Hartmann Date: Tue, 6 Mar 2001 11:51:58 +0000 Subject: [PATCH] Removed. git-svn-id: https://svn.dealii.org/trunk@4118 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/scripts/1d/lagrange | 89 ------- .../source/fe/scripts/1d/lagrange-cubic | 21 -- .../source/fe/scripts/1d/lagrange-quartic | 22 -- .../deal.II/source/fe/scripts/1d/postprocess | 29 -- .../source/fe/scripts/1d/restriction_dg | 100 ------- .../source/fe/scripts/1d/restriction_dg1 | 40 --- .../source/fe/scripts/1d/restriction_dg2 | 41 --- .../source/fe/scripts/1d/restriction_dg3 | 42 --- .../source/fe/scripts/1d/restriction_dg4 | 43 --- deal.II/deal.II/source/fe/scripts/2d/lagrange | 209 --------------- .../source/fe/scripts/2d/lagrange-cubic | 69 ----- .../source/fe/scripts/2d/lagrange-quadratic | 69 ----- .../source/fe/scripts/2d/lagrange-quartic | 83 ------ .../deal.II/source/fe/scripts/2d/postprocess | 42 --- .../source/fe/scripts/2d/restriction_dg | 121 --------- .../source/fe/scripts/2d/restriction_dg1 | 44 ---- .../source/fe/scripts/2d/restriction_dg2 | 40 --- .../source/fe/scripts/2d/restriction_dg3 | 58 ---- .../source/fe/scripts/2d/restriction_dg4 | 68 ----- .../fe/scripts/3d/computations_on_faces | 34 --- deal.II/deal.II/source/fe/scripts/3d/lagrange | 248 ------------------ .../source/fe/scripts/3d/lagrange-cubic | 161 ------------ .../source/fe/scripts/3d/lagrange-linear | 128 --------- .../source/fe/scripts/3d/lagrange-quadratic | 113 -------- .../source/fe/scripts/3d/lagrange-tools | 196 -------------- .../deal.II/source/fe/scripts/3d/postprocess | 60 ----- .../source/fe/scripts/3d/restriction_dg | 146 ----------- .../source/fe/scripts/3d/restriction_dg1 | 50 ---- .../source/fe/scripts/3d/restriction_dg2 | 59 ----- .../source/fe/scripts/3d/restriction_dg3 | 62 ----- 30 files changed, 2487 deletions(-) delete mode 100644 deal.II/deal.II/source/fe/scripts/1d/lagrange delete mode 100644 deal.II/deal.II/source/fe/scripts/1d/lagrange-cubic delete mode 100644 deal.II/deal.II/source/fe/scripts/1d/lagrange-quartic delete mode 100644 deal.II/deal.II/source/fe/scripts/1d/postprocess delete mode 100644 deal.II/deal.II/source/fe/scripts/1d/restriction_dg delete mode 100644 deal.II/deal.II/source/fe/scripts/1d/restriction_dg1 delete mode 100644 deal.II/deal.II/source/fe/scripts/1d/restriction_dg2 delete mode 100644 deal.II/deal.II/source/fe/scripts/1d/restriction_dg3 delete mode 100644 deal.II/deal.II/source/fe/scripts/1d/restriction_dg4 delete mode 100644 deal.II/deal.II/source/fe/scripts/2d/lagrange delete mode 100644 deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic delete mode 100644 deal.II/deal.II/source/fe/scripts/2d/lagrange-quadratic delete mode 100644 deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic delete mode 100644 deal.II/deal.II/source/fe/scripts/2d/postprocess delete mode 100644 deal.II/deal.II/source/fe/scripts/2d/restriction_dg delete mode 100644 deal.II/deal.II/source/fe/scripts/2d/restriction_dg1 delete mode 100644 deal.II/deal.II/source/fe/scripts/2d/restriction_dg2 delete mode 100644 deal.II/deal.II/source/fe/scripts/2d/restriction_dg3 delete mode 100644 deal.II/deal.II/source/fe/scripts/2d/restriction_dg4 delete mode 100644 deal.II/deal.II/source/fe/scripts/3d/computations_on_faces delete mode 100644 deal.II/deal.II/source/fe/scripts/3d/lagrange delete mode 100644 deal.II/deal.II/source/fe/scripts/3d/lagrange-cubic delete mode 100644 deal.II/deal.II/source/fe/scripts/3d/lagrange-linear delete mode 100644 deal.II/deal.II/source/fe/scripts/3d/lagrange-quadratic delete mode 100644 deal.II/deal.II/source/fe/scripts/3d/lagrange-tools delete mode 100644 deal.II/deal.II/source/fe/scripts/3d/postprocess delete mode 100644 deal.II/deal.II/source/fe/scripts/3d/restriction_dg delete mode 100644 deal.II/deal.II/source/fe/scripts/3d/restriction_dg1 delete mode 100644 deal.II/deal.II/source/fe/scripts/3d/restriction_dg2 delete mode 100644 deal.II/deal.II/source/fe/scripts/3d/restriction_dg3 diff --git a/deal.II/deal.II/source/fe/scripts/1d/lagrange b/deal.II/deal.II/source/fe/scripts/1d/lagrange deleted file mode 100644 index b57a028a09..0000000000 --- a/deal.II/deal.II/source/fe/scripts/1d/lagrange +++ /dev/null @@ -1,89 +0,0 @@ -# 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. 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); - grad_grad_phi_polynom := array(0..n_functions-1); - local_mass_matrix := array(0..n_functions-1, 0..n_functions-1); - - for i from 0 to n_functions-1 do - # note that the interp function wants vectors indexed from - # one and not from zero. - 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); - grad_phi_polynom[i] := diff(phi_polynom[i], xi); - grad_grad_phi_polynom[i] := diff(grad_phi_polynom[i], xi); - od; - - phi:= proc(i,x) 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; - - 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; - - - # to get the restriction (interpolation) matrices, evaluate - # the basis functions on the child cells at the global - # interpolation points - child_phi[0] := proc(i, point) - if ((point<0) or (point>1/2)) then - 0: - else - phi(i,2*point): - fi: - end: - child_phi[1] := proc(i, point) - if ((point<1/2) or (point>1)) then - 0: - else - phi(i,2*point-1): - fi: - end: - restriction := array(0..1,0..n_functions-1, 0..n_functions-1); - for child from 0 to 1 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]): - od: - od: - od: - - - for i from 0 to n_functions-1 do - for j from 0 to n_functions-1 do - local_mass_matrix[i,j] := int(phi_polynom[i] * phi_polynom[j] * h, - xi=0..1); - od; - od; - diff --git a/deal.II/deal.II/source/fe/scripts/1d/lagrange-cubic b/deal.II/deal.II/source/fe/scripts/1d/lagrange-cubic deleted file mode 100644 index a33cbd26f2..0000000000 --- a/deal.II/deal.II/source/fe/scripts/1d/lagrange-cubic +++ /dev/null @@ -1,21 +0,0 @@ - n_functions := 4; - - 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; - - - # do the real work - read "lagrange": - - - # write data to files - readlib(C); - C(phi_polynom, filename=cubic1d_shape_value); - C(grad_phi_polynom, filename=cubic1d_shape_grad); - C(grad_grad_phi_polynom, filename=cubic1d_shape_grad_grad); - C(prolongation, filename=cubic1d_prolongation); - C(restriction, filename=cubic1d_restriction); - C(local_mass_matrix, optimized, filename=cubic1d_massmatrix); diff --git a/deal.II/deal.II/source/fe/scripts/1d/lagrange-quartic b/deal.II/deal.II/source/fe/scripts/1d/lagrange-quartic deleted file mode 100644 index bf82c87237..0000000000 --- a/deal.II/deal.II/source/fe/scripts/1d/lagrange-quartic +++ /dev/null @@ -1,22 +0,0 @@ - n_functions := 5; - - 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; - - - # do the real work - read "lagrange": - - - # write data to files - readlib(C); - C(phi_polynom, filename=quartic1d_shape_value); - C(grad_phi_polynom, filename=quartic1d_shape_grad); - C(grad_grad_phi_polynom, filename=quartic1d_shape_grad_grad); - C(prolongation, filename=quartic1d_prolongation); - C(restriction, filename=quartic1d_restriction); - C(local_mass_matrix, optimized, filename=quartic1d_massmatrix); diff --git a/deal.II/deal.II/source/fe/scripts/1d/postprocess b/deal.II/deal.II/source/fe/scripts/1d/postprocess deleted file mode 100644 index 45c5371a60..0000000000 --- a/deal.II/deal.II/source/fe/scripts/1d/postprocess +++ /dev/null @@ -1,29 +0,0 @@ -# Use the following perl scripts to convert the output into a -# suitable format: -# -# $Id$ -# Wolfgang Bangerth, 1998 - - -# concatenate lines belonging together -perl -pi -e 's/([^;])\n/$1/g;' *1d_shape_value -perl -pi -e 's/([^;])\n/$1/g;' *1d_shape_grad -perl -pi -e 's/([^;])\n/$1/g;' *1d_shape_grad_grad - -# give the programs a structure -perl -pi -e 's/phi_polynom\[(\d)\] =/case $1: return/g;' *1d_shape_value -perl -pi -e 's/grad_phi_polynom\[(\d)\] = (.*);/case $1: return Point<1>($2);/g;' *1d_shape_grad -perl -pi -e 's/grad_grad_phi_polynom\[(\d)\] = (.*);/case $1: return_value[0][0] = $2;/g;' *1d_shape_grad_grad - -# use other indexing format for matrices -perl -pi -e 's/\[(\d+)\]\[(\d)\]/($1,$2)/g;' *1d_massmatrix -perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' *1d_prolongation -perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' *1d_restriction - -# give temporaries a data type -perl -pi -e 's/(t\d+) =/const double $1 =/g;' *1d_massmatrix - -# omit lines assigning zeroes to matrix elements, since zero is -# already set and to save compilation time -perl -pi -e 's/.*= 0.0;\n//g;' *1d_restriction -perl -pi -e 's/.*= 0.0;\n//g;' *1d_prolongation diff --git a/deal.II/deal.II/source/fe/scripts/1d/restriction_dg b/deal.II/deal.II/source/fe/scripts/1d/restriction_dg deleted file mode 100644 index 6bdcef8b25..0000000000 --- a/deal.II/deal.II/source/fe/scripts/1d/restriction_dg +++ /dev/null @@ -1,100 +0,0 @@ - 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 deleted file mode 100644 index 5e143fea43..0000000000 --- a/deal.II/deal.II/source/fe/scripts/1d/restriction_dg1 +++ /dev/null @@ -1,40 +0,0 @@ -# --------------------------------- 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 deleted file mode 100644 index 75639c3018..0000000000 --- a/deal.II/deal.II/source/fe/scripts/1d/restriction_dg2 +++ /dev/null @@ -1,41 +0,0 @@ -# --------------------------------- 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 deleted file mode 100644 index ba3dd1745e..0000000000 --- a/deal.II/deal.II/source/fe/scripts/1d/restriction_dg3 +++ /dev/null @@ -1,42 +0,0 @@ -# --------------------------------- 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 deleted file mode 100644 index b3b2c34817..0000000000 --- a/deal.II/deal.II/source/fe/scripts/1d/restriction_dg4 +++ /dev/null @@ -1,43 +0,0 @@ -# --------------------------------- 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/lagrange b/deal.II/deal.II/source/fe/scripts/2d/lagrange deleted file mode 100644 index a4352331f2..0000000000 --- a/deal.II/deal.II/source/fe/scripts/2d/lagrange +++ /dev/null @@ -1,209 +0,0 @@ -# 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. 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..1): - grad_grad_phi_polynom := array(0..n_functions-1,0..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): - - 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,1,0] := diff(phi_polynom[i], eta,xi): - grad_grad_phi_polynom[i,1,1] := diff(phi_polynom[i], eta,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. - # - # 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*. - 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): - Jacobian := linalg[matrix](2,2, [[diff(x_real,xi), diff(x_real,eta)], - [diff(y_real,xi), diff(y_real,eta)]]): - inverseJacobian := linalg[inverse](Jacobian): - detJ := linalg[det](Jacobian): - - grad_inverseJacobian := array(1..2, 1..2, 1..2): - for i from 1 to 2 do - for j from 1 to 2 do - for k from 1 to 2 do - if (i=1) then - grad_inverseJacobian[i,j,k] := diff (inverseJacobian[j,k], xi): - else - grad_inverseJacobian[i,j,k] := diff (inverseJacobian[j,k], eta): - fi: - od: - od: - od: - - 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 deleted file mode 100644 index 661117774d..0000000000 --- a/deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic +++ /dev/null @@ -1,69 +0,0 @@ -# --------------------------------- 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(grad_grad_phi_polynom, filename=cubic2d_shape_grad_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 deleted file mode 100644 index bcf0601c45..0000000000 --- a/deal.II/deal.II/source/fe/scripts/2d/lagrange-quadratic +++ /dev/null @@ -1,69 +0,0 @@ -# --------------------------------- 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(grad_grad_phi_polynom, filename=quadratic2d_shape_grad_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): - C(inverseJacobian, optimized, filename=quadratic2d_inverse_jacobian): - C(grad_inverseJacobian, optimized, - filename=quadratic2d_inverse_jacobian_grad): diff --git a/deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic b/deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic deleted file mode 100644 index 689ba343c9..0000000000 --- a/deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic +++ /dev/null @@ -1,83 +0,0 @@ -# --------------------------------- 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(grad_grad_phi_polynom, filename=quartic2d_shape_grad_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 deleted file mode 100644 index 84d6ad4c1c..0000000000 --- a/deal.II/deal.II/source/fe/scripts/2d/postprocess +++ /dev/null @@ -1,42 +0,0 @@ -# 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;' *2d_shape_value -perl -pi -e 's/xi/p(0)/g; s/zeta/p(2)/g; s/eta/p(1)/g;' *2d_shape_value - -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 - -perl -pi -e 's/^\s*t/const double t/g;' *2d_inverse_jacobian -perl -pi -e 's/x\[(\d)\]/vertices[$1](0)/g;' *2d_inverse_jacobian -perl -pi -e 's/y\[(\d)\]/vertices[$1](1)/g;' *2d_inverse_jacobian -perl -pi -e 's/inverseJacobian/jacobians[point]/g;' *2d_inverse_jacobian -perl -pi -e 's/\[(\d)\]\[(\d)\] =/($1,$2) =/g;' *2d_inverse_jacobian - -perl -pi -e 's/^\s*t/const double t/g;' *2d_inverse_jacobian_grad -perl -pi -e 's/x\[(\d)\]/vertices[$1](0)/g;' *2d_inverse_jacobian_grad -perl -pi -e 's/y\[(\d)\]/vertices[$1](1)/g;' *2d_inverse_jacobian_grad -perl -pi -e 's/inverseJacobian/jacobians_grad[point]/g;' *2d_inverse_jacobian_grad - -perl -pi -e 's/([^;])\n/$1/g;' *3d_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+)\]\[[12]\] = (.*);/$2);/g;' *2d_shape_grad -perl -pi -e 's/xi/p(0)/g; s/zeta/p(2)/g; s/eta/p(1)/g;' *2d_shape_grad - - -# concatenate all lines for each entry -perl -pi -e 's/([^;])\n/$1/g;' *2d_shape_grad_grad -# rename the variable -perl -pi -e 's/\s*grad_grad_phi_polynom/return_value/g;' *2d_shape_grad_grad -# insert 'case' and 'break' statements -perl -pi -e 's/(return_value\[(\d)\]\[0\]\[0\] = .*;)/break;\ncase $2:\n$1/g;' *2d_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;' *2d_shape_grad_grad diff --git a/deal.II/deal.II/source/fe/scripts/2d/restriction_dg b/deal.II/deal.II/source/fe/scripts/2d/restriction_dg deleted file mode 100644 index a789e31004..0000000000 --- a/deal.II/deal.II/source/fe/scripts/2d/restriction_dg +++ /dev/null @@ -1,121 +0,0 @@ - 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 deleted file mode 100644 index bf8c1da2ca..0000000000 --- a/deal.II/deal.II/source/fe/scripts/2d/restriction_dg1 +++ /dev/null @@ -1,44 +0,0 @@ -# --------------------------------- 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 deleted file mode 100644 index c093e67d9b..0000000000 --- a/deal.II/deal.II/source/fe/scripts/2d/restriction_dg2 +++ /dev/null @@ -1,40 +0,0 @@ -# --------------------------------- 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 deleted file mode 100644 index f4d60ec4ad..0000000000 --- a/deal.II/deal.II/source/fe/scripts/2d/restriction_dg3 +++ /dev/null @@ -1,58 +0,0 @@ -# --------------------------------- 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 deleted file mode 100644 index 3f4480053e..0000000000 --- a/deal.II/deal.II/source/fe/scripts/2d/restriction_dg4 +++ /dev/null @@ -1,68 +0,0 @@ -# --------------------------------- 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/computations_on_faces b/deal.II/deal.II/source/fe/scripts/3d/computations_on_faces deleted file mode 100644 index 0f80762cc9..0000000000 --- a/deal.II/deal.II/source/fe/scripts/3d/computations_on_faces +++ /dev/null @@ -1,34 +0,0 @@ - # 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 := sqrt (outward_vector[1]*outward_vector[1] + - outward_vector[2]*outward_vector[2] + - outward_vector[3]*outward_vector[3]): - normal_vector := map (p->p/detJ, outward_vector): - - measure := int (int (detJ, xi=0..1), eta=0..1): - - 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 deleted file mode 100644 index 7dc2549e7b..0000000000 --- a/deal.II/deal.II/source/fe/scripts/3d/lagrange +++ /dev/null @@ -1,248 +0,0 @@ -# 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-cubic b/deal.II/deal.II/source/fe/scripts/3d/lagrange-cubic deleted file mode 100644 index 1716d22041..0000000000 --- a/deal.II/deal.II/source/fe/scripts/3d/lagrange-cubic +++ /dev/null @@ -1,161 +0,0 @@ -# --------------------------------- 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 := 64: - n_face_functions := 16: - n_constraints := 45: - - 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: - 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, 2, support_points): - support_points_fill_quads (32, 2, support_points): - support_points_fill_hex (56, 2, support_points): - - 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/3,0]: - face_support_points[5] := [2/3,0]: - face_support_points[6] := [1,1/3]: - face_support_points[7] := [1,2/3]: - face_support_points[8] := [1/3,1]: - face_support_points[9] := [2/3,1]: - face_support_points[10] := [0,1/3]: - face_support_points[11] := [0,2/3]: - face_support_points[12] := [1/3,1/3]: - face_support_points[13] := [1/3,2/3]: - face_support_points[14] := [2/3,1/3]: - face_support_points[15] := [2/3,2/3]: - - # 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] := 12: - constrained_face_function[9] := 13: - constrained_face_function[10] := 14: - constrained_face_function[11] := 15: - constrained_face_function[12] := 32: - constrained_face_function[13] := 33: - constrained_face_function[14] := 34: - constrained_face_function[15] := 35: - - 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/6]): # lines from center to boundary - constrained_face_support_points[6] := array(0..1, [1/2,2/6]): - constrained_face_support_points[7] := array(0..1, [4/6,1/2]): - constrained_face_support_points[8] := array(0..1, [5/6,1/2]): - constrained_face_support_points[9] := array(0..1, [1/2,4/6]): - constrained_face_support_points[10] := array(0..1, [1/2,5/6]): - constrained_face_support_points[11] := array(0..1, [1/6,1/2]): - constrained_face_support_points[12] := array(0..1, [2/6,1/2]): - constrained_face_support_points[13] := array(0..1, [1/6,0]): # children of bounding lines - constrained_face_support_points[14] := array(0..1, [2/6,0]): - constrained_face_support_points[15] := array(0..1, [4/6,0]): - constrained_face_support_points[16] := array(0..1, [5/6,0]): - - constrained_face_support_points[17] := array(0..1, [1,1/6]): - constrained_face_support_points[18] := array(0..1, [1,2/6]): - constrained_face_support_points[19] := array(0..1, [1,4/6]): - constrained_face_support_points[20] := array(0..1, [1,5/6]): - - constrained_face_support_points[21] := array(0..1, [1/6,1]): - constrained_face_support_points[22] := array(0..1, [2/6,1]): - constrained_face_support_points[23] := array(0..1, [4/6,1]): - constrained_face_support_points[24] := array(0..1, [5/6,1]): - - constrained_face_support_points[25] := array(0..1, [0,1/6]): - constrained_face_support_points[26] := array(0..1, [0,2/6]): - constrained_face_support_points[27] := array(0..1, [0,4/6]): - constrained_face_support_points[28] := array(0..1, [0,5/6]): - - constrained_face_support_points[29] := array(0..1, [1/6,1/6]): # child quads - constrained_face_support_points[30] := array(0..1, [2/6,1/6]): - constrained_face_support_points[31] := array(0..1, [1/6,2/6]): - constrained_face_support_points[32] := array(0..1, [2/6,2/6]): - - constrained_face_support_points[33] := array(0..1, [4/6,1/6]): - constrained_face_support_points[34] := array(0..1, [5/6,1/6]): - constrained_face_support_points[35] := array(0..1, [4/6,2/6]): - constrained_face_support_points[36] := array(0..1, [5/6,2/6]): - - constrained_face_support_points[37] := array(0..1, [4/6,4/6]): - constrained_face_support_points[38] := array(0..1, [5/6,4/6]): - constrained_face_support_points[39] := array(0..1, [4/6,5/6]): - constrained_face_support_points[40] := array(0..1, [5/6,5/6]): - - constrained_face_support_points[41] := array(0..1, [1/6,4/6]): - constrained_face_support_points[42] := array(0..1, [2/6,4/6]): - constrained_face_support_points[43] := array(0..1, [1/6,5/6]): - constrained_face_support_points[44] := array(0..1, [2/6,5/6]): - - - # do the real work - read "lagrange": - - - - # write data to files - print ("writing data to files"): - readlib(C): - C(phi_polynom, filename=cubic3d_shape_value): - C(grad_phi_polynom, filename=cubic3d_shape_grad): - C(grad_grad_phi_polynom, filename=cubic3d_shape_grad_grad): - C(prolongation, filename=cubic3d_prolongation): - C(restriction, filename=cubic3d_restriction): - C(interface_constraints, filename=cubic3d_constraints): - C(real_space_points, optimized, filename=cubic3d_real_points): - - writeto (cubic3d_unit_support_points): - print (support_points): - - \ 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 deleted file mode 100644 index afcc425b59..0000000000 --- a/deal.II/deal.II/source/fe/scripts/3d/lagrange-linear +++ /dev/null @@ -1,128 +0,0 @@ -# --------------------------------- 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): - # the following two files get much smaller and faster when processed using 'optimized' - C(inverseJacobian, optimized, filename=linear3d_inverse_jacobian): - C(grad_inverseJacobian, optimized, 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 deleted file mode 100644 index 69dfc0adda..0000000000 --- a/deal.II/deal.II/source/fe/scripts/3d/lagrange-quadratic +++ /dev/null @@ -1,113 +0,0 @@ -# --------------------------------- 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 deleted file mode 100644 index 9a27c773b2..0000000000 --- a/deal.II/deal.II/source/fe/scripts/3d/lagrange-tools +++ /dev/null @@ -1,196 +0,0 @@ -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 - support_points[next_index] - := 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] - := 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] - := 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] - := 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] - := 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] - := 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] - := 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] - := 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] - := 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] - := 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] - := 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] - := array (1..3, [0, i*increment, 1]): - next_index := next_index+1 - od: -end: - - - - -support_points_fill_quads := - proc (starting_index, dofs_per_direction, support_points) - - local next_index, increment, i,j: - - next_index := starting_index: - increment := 1/(dofs_per_direction+1): - - # face 0 - for i from 1 to dofs_per_direction do - for j from 1 to dofs_per_direction do - support_points[next_index] - := array (1..3, [j*increment, 0, i*increment]): - next_index := next_index+1: - od: - od: - - # face 1 - for i from 1 to dofs_per_direction do - for j from 1 to dofs_per_direction do - support_points[next_index] - := array (1..3, [j*increment, 1, i*increment]): - next_index := next_index+1: - od: - od: - - # face 2 - for i from 1 to dofs_per_direction do - for j from 1 to dofs_per_direction do - support_points[next_index] - := array (1..3, [j*increment, i*increment, 0]): - next_index := next_index+1: - od: - od: - - # face 3 - for i from 1 to dofs_per_direction do - for j from 1 to dofs_per_direction do - support_points[next_index] - := array (1..3, [1, j*increment, i*increment]): - next_index := next_index+1: - od: - od: - - - - # face 4 - for i from 1 to dofs_per_direction do - for j from 1 to dofs_per_direction do - support_points[next_index] - := array (1..3, [j*increment, i*increment, 1]): - next_index := next_index+1: - od: - od: - - # face 5 - for i from 1 to dofs_per_direction do - for j from 1 to dofs_per_direction do - support_points[next_index] - := array (1..3, [0, j*increment, i*increment]): - next_index := next_index+1: - od: - od: -end: - - - - -support_points_fill_hex := proc (starting_index, dofs_per_direction, support_points) - - local next_index, increment, i, j, k: - - next_index := starting_index: - increment := 1/(dofs_per_direction+1): - - for i from 1 to dofs_per_direction do - for j from 1 to dofs_per_direction do - for k from 1 to dofs_per_direction do - support_points[next_index] - := array (1..3, [k*increment, - j*increment, - i*increment]): - next_index := next_index + 1: - od: - od: - od: -end: diff --git a/deal.II/deal.II/source/fe/scripts/3d/postprocess b/deal.II/deal.II/source/fe/scripts/3d/postprocess deleted file mode 100644 index 1f26bb73e1..0000000000 --- a/deal.II/deal.II/source/fe/scripts/3d/postprocess +++ /dev/null @@ -1,60 +0,0 @@ -# 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/([^;])\n/$1/g;' *3d_inverse_jacobian -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/z\[(\d)\]/vertices[$1](2)/g;' *3d_inverse_jacobian -perl -pi -e 's/^s*inverseJacobian/ jacobians[point]/g;' *3d_inverse_jacobian -perl -pi -e 's/\[(\d)\]\[(\d)\] =/($1,$2) =/g;' *3d_inverse_jacobian - -perl -pi -e 's/([^;])\n/$1/g;' *3d_inverse_jacobian_grad -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/z\[(\d)\]/vertices[$1](2)/g;' *3d_inverse_jacobian_grad -perl -pi -e 's/^\s*grad_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 diff --git a/deal.II/deal.II/source/fe/scripts/3d/restriction_dg b/deal.II/deal.II/source/fe/scripts/3d/restriction_dg deleted file mode 100644 index 9bcc76ec55..0000000000 --- a/deal.II/deal.II/source/fe/scripts/3d/restriction_dg +++ /dev/null @@ -1,146 +0,0 @@ - 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 deleted file mode 100644 index 05dae43b0c..0000000000 --- a/deal.II/deal.II/source/fe/scripts/3d/restriction_dg1 +++ /dev/null @@ -1,50 +0,0 @@ -# --------------------------------- 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 deleted file mode 100644 index 14dcd93dcb..0000000000 --- a/deal.II/deal.II/source/fe/scripts/3d/restriction_dg2 +++ /dev/null @@ -1,59 +0,0 @@ -# --------------------------------- 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 deleted file mode 100644 index a315aa5422..0000000000 --- a/deal.II/deal.II/source/fe/scripts/3d/restriction_dg3 +++ /dev/null @@ -1,62 +0,0 @@ -# --------------------------------- 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); - - - - - - - - - - - -- 2.39.5