]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Extend scripts slowly for the calculation of second derivatives.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 4 Nov 1998 19:05:41 +0000 (19:05 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 4 Nov 1998 19:05:41 +0000 (19:05 +0000)
git-svn-id: https://svn.dealii.org/trunk@632 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/fe/scripts/2d/lagrange
deal.II/deal.II/source/fe/scripts/2d/lagrange-cubic
deal.II/deal.II/source/fe/scripts/2d/lagrange-quadratic
deal.II/deal.II/source/fe/scripts/2d/lagrange-quartic
deal.II/deal.II/source/fe/scripts/2d/postprocess

index b66816fe0ec4b23b5f2b618822a8a682da1db99e..a4352331f2236c8a06c6da3dc6b59c9ca8f4db43 100644 (file)
@@ -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,
     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:
   # 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[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
index 431b6781c7760ad6601d7605f19f1cc817c767e4..661117774de37d978b5f3f90c5f0c1e80a075fb9 100644 (file)
   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):
 
index deb54046223c07d072c5ac4451604155149e6cf6..bcf0601c4501a62c4346899ee1f6f399c6b4522f 100644 (file)
 
 
   # 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):
index d6ced17de01623edd7d4408345b2d89d279ac4fc..689ba343c9916e74eee6d10cc5bb9b2a8f7e3272 100644 (file)
 
 
   # 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):
 
index 9b06f5cbcb94077b781a7f7c90fc462271a1bc20..b56e78696382a28cc3f4847d6f353e624b0c5f43 100644 (file)
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.