+++ /dev/null
-# 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;
-
+++ /dev/null
- 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);
+++ /dev/null
- 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);
+++ /dev/null
-# 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
+++ /dev/null
- 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:
-
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-# --------------------------------- 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);
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-# --------------------------------- 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);
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-# --------------------------------- 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);
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-# --------------------------------- 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);
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-# 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:
+++ /dev/null
-# --------------------------------- 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):
-
+++ /dev/null
-# --------------------------------- 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):
+++ /dev/null
-# --------------------------------- 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):
-
+++ /dev/null
-# 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
+++ /dev/null
- 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:
-
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-# --------------------------------- 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);
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-# --------------------------------- 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);
-
+++ /dev/null
-# --------------------------------- 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);
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-# --------------------------------- 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);
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
- # 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
+++ /dev/null
-# 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
+++ /dev/null
-# --------------------------------- 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
+++ /dev/null
-# --------------------------------- 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):
-
-
+++ /dev/null
-# --------------------------------- 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
+++ /dev/null
-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:
+++ /dev/null
-# 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
+++ /dev/null
- 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:
-
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-# --------------------------------- 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);
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-# --------------------------------- 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);
-
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-# --------------------------------- 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);
-
-
-
-
-
-
-
-
-
-
-