From: wolf Date: Wed, 4 Nov 1998 19:05:41 +0000 (+0000) Subject: Extend scripts slowly for the calculation of second derivatives. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0a56a6f350dd316da06e30fe89949cbb38789f11;p=dealii-svn.git Extend scripts slowly for the calculation of second derivatives. git-svn-id: https://svn.dealii.org/trunk@632 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/scripts/2d/lagrange b/deal.II/deal.II/source/fe/scripts/2d/lagrange index b66816fe0e..a4352331f2 100644 --- a/deal.II/deal.II/source/fe/scripts/2d/lagrange +++ b/deal.II/deal.II/source/fe/scripts/2d/lagrange @@ -1,12 +1,14 @@ # 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. +# 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, @@ -37,6 +39,11 @@ 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: @@ -122,6 +129,14 @@ # 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): @@ -130,7 +145,24 @@ 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): + 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 diff --git a/deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic b/deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic index 431b6781c7..661117774d 100644 --- a/deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic +++ b/deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic @@ -52,17 +52,18 @@ constrained_face_support_points[4] := 5/6: # do the real work - read "lagrange" + 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): + 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 index deb5404622..bcf0601c45 100644 --- a/deal.II/deal.II/source/fe/scripts/2d/lagrange-quadratic +++ b/deal.II/deal.II/source/fe/scripts/2d/lagrange-quadratic @@ -50,17 +50,20 @@ # do the real work - read "lagrange" + 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): - + 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 index d6ced17de0..689ba343c9 100644 --- a/deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic +++ b/deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic @@ -66,17 +66,18 @@ # do the real work - read "lagrange" + 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): + 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 index 9b06f5cbcb..b56e786963 100644 --- a/deal.II/deal.II/source/fe/scripts/2d/postprocess +++ b/deal.II/deal.II/source/fe/scripts/2d/postprocess @@ -4,14 +4,19 @@ # $Id$ # Wolfgang Bangerth, 1998 -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 +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 +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 \ No newline at end of file