* Return a readonly reference to the
* matrix which describes the transfer of a
* child with the given number to the
- * mother cell.
+ * mother cell. See the #restriction# array
+ * for more information.
*/
const dFMatrix & restrict (const unsigned int child) const;
* are for the refined cell's degrees of
* freedom.
*
+ * In essence, using the matrices from the
+ * children to the mother cell amounts to
+ * computing the interpolation of the
+ * function on the refined to the coarse
+ * mesh. To get the vector of nodal values
+ * of the interpolant on the mother cell,
+ * you have to multiply the nodal value
+ * vectors of each of the child cell with
+ * the respective restriction matrix and
+ * clobber these contributions together.
+ * However, you must take care not to
+ * #add# these together, since nodes which
+ * belong to more than one child would then
+ * be counted more than once; rather, you
+ * have to overwrite the nonzero
+ * contributions of each child into the
+ * nodal value vector of the mother cell.
+ *
+ * To compute the interpolation of a
+ * finite element field to a cell, you
+ * may use the #interpolation# function.
+ * See there for more information.
+ *
* Upon assembling the transfer matrix
* between cells using this matrix array,
* zero elements in the restriction
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:
+
+
# these are the basis functions differentiated with respect to
# xi and eta. we need them for the computation of the jacobi
# matrix, since we can't just differentiate a function.
print ("writing data to files"):
readlib(C):
C(prolongation, filename=prolongation_2d):
+ C(restriction, filename=restriction_2d):
C(array(1..2, [x[4], y[4]]), optimized, filename=crosspoint_2d):
C(mass_matrix, optimized, filename=massmatrix_2d):
prolongation[3](3,3) = 1.0;
prolongation[3](4,3) = 1.0/2.0;
prolongation[3](4,4) = 1.0/2.0;
+
+ restriction[0](0,0) = 1.0;
+ restriction[0](4,2) = 1.0;
+ restriction[1](1,1) = 1.0;
+ restriction[1](4,3) = 1.0;
+ restriction[2](2,2) = 1.0;
+ restriction[2](4,0) = 1.0;
+ restriction[3](3,3) = 1.0;
+ restriction[3](4,1) = 1.0;
};
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,
C(phi_polynom, filename=shape_value_1d);
C(grad_phi_polynom, filename=shape_grad_1d);
C(prolongation, filename=prolongation_1d);
+ C(restriction, filename=restriction_1d);
C(local_mass_matrix, optimized, filename=massmatrix_1d);
-----------------------------------------------------------------------
perl -pi -e 's/grad_phi_polynom\[(\d)\] = (.*);/case $1: return Point<1>($2);/g;' shape_grad_1d
perl -pi -e 's/\[(\d+)\]\[(\d)\]/($1,$2)/g;' massmatrix_1d
perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' prolongation_1d
+ perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' restriction_1d
perl -pi -e 's/(t\d+) =/const double $1 =/g;' massmatrix_1d
*/
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
C(phi_polynom, filename=shape_value_2d):
C(grad_phi_polynom, filename=shape_grad_2d):
C(prolongation, filename=prolongation_2d):
+ C(restriction, filename=restriction_2d):
C(local_mass_matrix, optimized, filename=massmatrix_2d):
C(interface_constraints, filename=constraints_2d):
C(real_points, optimized, filename=real_points_2d);
perl -pi -e 's/(t\d+) =/const double $1 =/g;' massmatrix_2d
perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' prolongation_2d
perl -pi -e 's/.*= 0.0;\n//g;' prolongation_2d
+ perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' restriction_2d
+ perl -pi -e 's/.*= 0.0;\n//g;' restriction_2d
perl -pi -e 's/\[(\d+)\]\[(\d+)\]/($1,$2)/g;' constraints_2d
*/
prolongation[1](3,1) = 5.0/16.0;
prolongation[1](3,2) = -5.0/16.0;
prolongation[1](3,3) = 15.0/16.0;
+
+ restriction[0](0,0) = 1.0;
+ restriction[0](0,1) = 0.0;
+ restriction[0](0,2) = 0.0;
+ restriction[0](0,3) = 0.0;
+ restriction[0](1,0) = 0.0;
+ restriction[0](1,1) = 0.0;
+ restriction[0](1,2) = 0.0;
+ restriction[0](1,3) = 0.0;
+ restriction[0](2,0) = 0.0;
+ restriction[0](2,1) = 0.0;
+ restriction[0](2,2) = 0.0;
+ restriction[0](2,3) = 1.0;
+ restriction[0](3,0) = 0.0;
+ restriction[0](3,1) = 0.0;
+ restriction[0](3,2) = 0.0;
+ restriction[0](3,3) = 0.0;
+ restriction[1](0,0) = 0.0;
+ restriction[1](0,1) = 0.0;
+ restriction[1](0,2) = 0.0;
+ restriction[1](0,3) = 0.0;
+ restriction[1](1,0) = 0.0;
+ restriction[1](1,1) = 1.0;
+ restriction[1](1,2) = 0.0;
+ restriction[1](1,3) = 0.0;
+ restriction[1](2,0) = 0.0;
+ restriction[1](2,1) = 0.0;
+ restriction[1](2,2) = 0.0;
+ restriction[1](2,3) = 0.0;
+ restriction[1](3,0) = 0.0;
+ restriction[1](3,1) = 0.0;
+ restriction[1](3,2) = 1.0;
+ restriction[1](3,3) = 0.0;
};
prolongation[3](15,13) = 25.0/256.0;
prolongation[3](15,14) = -75.0/256.0;
prolongation[3](15,15) = 225.0/256.0;
+
+ restriction[0](0,0) = 1.0;
+ restriction[0](4,5) = 1.0;
+ restriction[0](10,11) = 1.0;
+ restriction[0](12,14) = 1.0;
+ restriction[1](1,1) = 1.0;
+ restriction[1](5,4) = 1.0;
+ restriction[1](6,7) = 1.0;
+ restriction[1](13,15) = 1.0;
+ restriction[2](2,2) = 1.0;
+ restriction[2](7,6) = 1.0;
+ restriction[2](9,8) = 1.0;
+ restriction[2](14,12) = 1.0;
+ restriction[3](3,3) = 1.0;
+ restriction[3](8,9) = 1.0;
+ restriction[3](11,10) = 1.0;
+ restriction[3](15,13) = 1.0;
};
// we do not add up the contributions but
// set them right into the matrices!
restriction[0](0,0) = 1.0;
- restriction[0](0,1) = 1./2.;
- restriction[0](1,1) = 1./2.;
-
- restriction[1](0,0) = 1./2.;
- restriction[1](1,0) = 1./2.;
restriction[1](1,1) = 1.0;
prolongation[0](0,0) = 1.0;
interface_constraints(0,1) = 1./2.;
restriction[0](0,0) = 1.0;
- restriction[0](0,1) = 1./2.;
- restriction[0](1,1) = 1./2.;
- restriction[0](0,3) = 1./2.;
- restriction[0](3,3) = 1./2.;
- restriction[0](0,2) = 1./4.;
- restriction[0](1,2) = 1./4.;
- restriction[0](2,2) = 1./4.;
- restriction[0](3,2) = 1./4.;
-
restriction[1](1,1) = 1.0;
- restriction[1](0,0) = 1./2.;
- restriction[1](1,0) = 1./2.;
- restriction[1](1,2) = 1./2.;
- restriction[1](2,2) = 1./2.;
- restriction[1](0,3) = 1./4.;
- restriction[1](1,3) = 1./4.;
- restriction[1](2,3) = 1./4.;
- restriction[1](3,3) = 1./4.;
-
restriction[2](2,2) = 1.0;
- restriction[2](2,1) = 1./2.;
- restriction[2](1,1) = 1./2.;
- restriction[2](2,3) = 1./2.;
- restriction[2](3,3) = 1./2.;
- restriction[2](0,0) = 1./4.;
- restriction[2](1,0) = 1./4.;
- restriction[2](2,0) = 1./4.;
- restriction[2](3,0) = 1./4.;
-
restriction[3](3,3) = 1.0;
- restriction[3](0,0) = 1./2.;
- restriction[3](3,0) = 1./2.;
- restriction[3](2,2) = 1./2.;
- restriction[3](3,2) = 1./2.;
- restriction[3](0,1) = 1./4.;
- restriction[3](1,1) = 1./4.;
- restriction[3](2,1) = 1./4.;
- restriction[3](3,1) = 1./4.;
prolongation[0](0,0) = 1.0;
prolongation[0](1,0) = 1./2.;
+/*-----------------------------------------------------------------
+ * For the 2D stuff, you may use the script below. However, 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.
+ 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:
+
+ phi_polynom := array(0..n_functions-1):
+ grad_phi_polynom := array(0..n_functions-1,0..1):
+ local_mass_matrix := array(0..n_functions-1, 0..n_functions-1):
+ prolongation := array(0..3,0..n_functions-1, 0..n_functions-1):
+ interface_constraints := array(0..2*(n_face_functions-2)+1-1,
+ 0..n_face_functions-1):
+ real_points := array(0..n_functions-1, 0..1);
+
+ print ("Computing basis functions"):
+ for i from 0 to n_functions-1 do
+ print (i):
+ values := array(1..n_functions):
+ for j from 1 to n_functions do
+ values[j] := 0:
+ od:
+ values[i+1] := 1:
+
+ equation_system := {}:
+ for j from 0 to n_functions-1 do
+ poly := subs(xi=support_points[j][1],
+ eta=support_points[j][2],
+ trial_function):
+ if (i=j) then
+ equation_system := equation_system union {poly = 1}:
+ else
+ equation_system := equation_system union {poly = 0}:
+ fi:
+ od:
+
+ phi_polynom[i] := subs(solve(equation_system), trial_function):
+ grad_phi_polynom[i,0] := diff(phi_polynom[i], xi):
+ grad_phi_polynom[i,1] := diff(phi_polynom[i], eta):
+ od:
+
+ phi:= proc(i,x,y) subs(xi=x, eta=y, phi_polynom[i]): end:
+
+
+ #points on children: let them be indexed one-based, as are
+ #the support_points
+ points[0] := array(0..n_functions-1, 1..2):
+ points[1] := array(0..n_functions-1, 1..2):
+ points[2] := array(0..n_functions-1, 1..2):
+ points[3] := array(0..n_functions-1, 1..2):
+ for i from 0 to n_functions-1 do
+ points[0][i,1] := support_points[i][1]/2:
+ points[0][i,2] := support_points[i][2]/2:
+
+ points[1][i,1] := support_points[i][1]/2+1/2:
+ points[1][i,2] := support_points[i][2]/2:
+
+ points[2][i,1] := support_points[i][1]/2+1/2:
+ points[2][i,2] := support_points[i][2]/2+1/2:
+
+ points[3][i,1] := support_points[i][1]/2:
+ points[3][i,2] := support_points[i][2]/2+1/2:
+ od:
+
+ print ("Computing prolongation matrices"):
+ for i from 0 to 3 do
+ print ("child", i):
+ for j from 0 to n_functions-1 do
+ for k from 0 to n_functions-1 do
+ prolongation[i,j,k] := phi(k, points[i][j,1], points[i][j,2]):
+ od:
+ od:
+ od:
+
+ print ("Computing restriction matrices"):
+ # to get the restriction (interpolation) matrices, evaluate
+ # the basis functions on the child cells at the global
+ # interpolation points
+ child_phi[0] := proc(i, x, y)
+ if ((x>1/2) or (y>1/2)) then
+ 0:
+ else
+ phi(i,2*x,2*y):
+ fi:
+ end:
+ child_phi[1] := proc(i, x, y)
+ if ((x<1/2) or (y>1/2)) then
+ 0:
+ else
+ phi(i,2*x-1,2*y):
+ fi:
+ end:
+ child_phi[2] := proc(i, x, y)
+ if ((x<1/2) or (y<1/2)) then
+ 0:
+ else
+ phi(i,2*x-1,2*y-1):
+ fi:
+ end:
+ child_phi[3] := proc(i, x, y)
+ if ((x>1/2) or (y<1/2)) then
+ 0:
+ else
+ phi(i,2*x,2*y-1):
+ fi:
+ end:
+ restriction := array(0..3,0..n_functions-1, 0..n_functions-1);
+ for child from 0 to 3 do
+ for j from 0 to n_functions-1 do
+ for k from 0 to n_functions-1 do
+ restriction[child,j,k] := child_phi[child](k,
+ support_points[j][1],
+ support_points[j][2]):
+ od:
+ od:
+ od:
+
+
+ print ("Computing local mass matrix"):
+ # tphi are the basis functions of the linear element. These functions
+ # are used for the computation of the subparametric transformation from
+ # unit cell to real cell.
+ # x and y are arrays holding the x- and y-values of the four vertices
+ # of this cell in real space.
+ x := array(0..3);
+ y := array(0..3);
+ tphi[0] := (1-xi)*(1-eta):
+ tphi[1] := xi*(1-eta):
+ tphi[2] := xi*eta:
+ tphi[3] := (1-xi)*eta:
+ x_real := sum(x[s]*tphi[s], s=0..3):
+ y_real := sum(y[s]*tphi[s], s=0..3):
+ detJ := diff(x_real,xi)*diff(y_real,eta) - diff(x_real,eta)*diff(y_real,xi):
+ for i from 0 to n_functions-1 do
+ print ("line", i):
+ for j from 0 to n_functions-1 do
+ local_mass_matrix[i,j] := int(int(phi_polynom[i] * phi_polynom[j] * detJ,
+ xi=0..1), eta=0..1):
+ od:
+ od:
+
+ print ("computing support points in real space"):
+ for i from 0 to n_functions-1 do
+ real_points[i,0] := subs(xi=support_points[i][1],
+ eta=support_points[i][2], x_real);
+ real_points[i,1] := subs(xi=support_points[i][1],
+ eta=support_points[i][2], y_real);
+ od:
+
+ print ("computing interface constraint matrices"):
+ # compute the interface constraint matrices. these are the values of the
+ # basis functions on the coarse cell's face at the points of the child
+ # cell's basis functions on the child faces
+ face_phi_polynom := array(0..n_face_functions-1):
+ for i from 0 to n_face_functions-1 do
+ # note that the interp function wants vectors indexed from
+ # one and not from zero.
+ values := array(1..n_face_functions):
+ for j from 1 to n_face_functions do
+ values[j] := 0:
+ od:
+ values[i+1] := 1:
+
+ shifted_face_support_points := array (1..n_face_functions):
+ for j from 1 to n_face_functions do
+ shifted_face_support_points[j] := face_support_points[j-1]:
+ od:
+
+ face_phi_polynom[i] := interp (shifted_face_support_points, values, xi):
+ od:
+
+ for i from 0 to 2*(n_face_functions-2)+1-1 do
+ for j from 0 to n_face_functions-1 do
+ interface_constraints[i,j] := subs(xi=constrained_face_support_points[i],
+ face_phi_polynom[j]);
+ od:
+ od:
+
+
+ print ("writing data to files"):
+ readlib(C):
+ C(phi_polynom, filename=shape_value_2d):
+ C(grad_phi_polynom, filename=shape_grad_2d):
+ C(prolongation, filename=prolongation_2d):
+ C(restriction, filename=restriction_2d):
+ C(local_mass_matrix, optimized, filename=massmatrix_2d):
+ C(interface_constraints, filename=constraints_2d):
+ C(real_points, optimized, filename=real_points_2d);
+
+ ---------------------------------------------------------------
+
+ Use the following perl scripts to convert the output into a
+ suitable format.
+
+ perl -pi -e 's/phi_polynom\[(\d+)\] =/case $1: return/g;' shape_value_2d
+ perl -pi -e 's/([^;])\n/$1/g;' shape_grad_2d
+ perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[0\] = (.*);/case $1: return Point<2>($2,/g;' shape_grad_2d
+ perl -pi -e 's/grad_phi_polynom\[(\d+)\]\[1\] = (.*);/$2);/g;' shape_grad_2d
+ perl -pi -e 's/\[(\d+)\]\[(\d+)\]/($1,$2)/g;' massmatrix_2d
+ perl -pi -e 's/(t\d+) =/const double $1 =/g;' massmatrix_2d
+ perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' prolongation_2d
+ perl -pi -e 's/.*= 0.0;\n//g;' prolongation_2d
+ perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' restriction_2d
+ perl -pi -e 's/.*= 0.0;\n//g;' restriction_2d
+ perl -pi -e 's/\[(\d+)\]\[(\d+)\]/($1,$2)/g;' constraints_2d
+*/
+
+
+
#if deal_II_dimension == 1
points[1] := array(0..2, [1/2, 1, 3/4]);
prolongation := array(0..1,0..2, 0..2);
+ restriction := array(0..1,0..2, 0..2);
for i from 0 to 1 do
for j from 0 to 2 do
od;
od;
+
+
+ # to get the restriction (interpolation) matrices, evaluate
+ # the basis functions on the child cells at the global
+ # interpolation points
+
+ global_points := array(0..2, [0,1,1/2]):
+ 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:
+
+ for child from 0 to 1 do
+ for j from 0 to 2 do
+ for k from 0 to 2 do
+ restriction[child,j,k] := child_phi[child](k, global_points[j]):
+ od:
+ od:
+ od:
+
readlib(C);
C(prolongation);
+ C(restriction);
*/
prolongation[0](0,0) = 1.0;
prolongation[1](2,0) = -1.0/8.0;
prolongation[1](2,1) = 3.0/8.0;
prolongation[1](2,2) = 3.0/4.0;
+
+ restriction[0](0,0)= 1.0;
+ restriction[0](2,1)= 1.0;
+ restriction[1](1,1)= 1.0;
+ restriction[1](2,0)= 1.0;
};
interface_constraints(2,2) = 3./4.;
/*
- Get the prolongation matrices by the following little maple script:
+ Get the prolongation and restriction matrices by the following little maple script:
phi[0] := proc(xi,eta) (1-xi)*( 2*xi-1) * (1-eta)*( 2*eta-1); end;
phi[1] := proc(xi,eta) xi *(-2*xi+1) * (1-eta)*( 2*eta-1); end;
prolongation[3](8,6) = 9.0/32.0;
prolongation[3](8,7) = 9.0/32.0;
prolongation[3](8,8) = 9.0/16.0;
+
+ restriction[0](0,0) = 1.0;
+ restriction[0](4,1) = 1.0;
+ restriction[0](7,3) = 1.0;
+ restriction[0](8,2) = 1.0;
+ restriction[1](1,1) = 1.0;
+ restriction[1](4,0) = 1.0;
+ restriction[1](5,2) = 1.0;
+ restriction[1](8,3) = 1.0;
+ restriction[2](2,2) = 1.0;
+ restriction[2](5,1) = 1.0;
+ restriction[2](6,3) = 1.0;
+ restriction[2](8,0) = 1.0;
+ restriction[3](3,3) = 1.0;
+ restriction[3](6,2) = 1.0;
+ restriction[3](7,0) = 1.0;
+ restriction[3](8,1) = 1.0;
};
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:
+
+
readlib(C);
C(phi_polynom, filename=shape_value_1d);
C(grad_phi_polynom, filename=shape_grad_1d);
C(prolongation, filename=prolongation_1d);
+ C(restriction, filename=restriction_1d);
C(local_mass_matrix, optimized, filename=massmatrix_1d);
-----------------------------------------------------------------------
perl -pi -e 's/\[(\d+)\]\[(\d)\]/($1,$2)/g;' massmatrix_1d
perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' prolongation_1d
perl -pi -e 's/.*= 0.0;\n//g;' prolongation_1d
+ perl -pi -e 's/\[(\d+)\]\[(\d)\]\[(\d)\]/[$1]($2,$3)/g;' restriction_1d
+ perl -pi -e 's/.*= 0.0;\n//g;' restriction_1d
perl -pi -e 's/(t\d+) =/const double $1 =/g;' massmatrix_1d
*/
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
C(phi_polynom, filename=shape_value_2d):
C(grad_phi_polynom, filename=shape_grad_2d):
C(prolongation, filename=prolongation_2d):
+ C(restriction, filename=restriction_2d);
C(local_mass_matrix, optimized, filename=massmatrix_2d):
C(interface_constraints, filename=constraints_2d):
C(real_points, optimized, filename=real_points_2d);
perl -pi -e 's/(t\d+) =/const double $1 =/g;' massmatrix_2d
perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' prolongation_2d
perl -pi -e 's/.*= 0.0;\n//g;' prolongation_2d
+ perl -pi -e 's/\[(\d+)\]\[(\d+)\]\[(\d+)\]/[$1]($2,$3)/g;' restriction_2d
+ perl -pi -e 's/.*= 0.0;\n//g;' restriction_2d
perl -pi -e 's/\[(\d+)\]\[(\d+)\]/($1,$2)/g;' constraints_2d
*/
prolongation[1](4,2) = 7.0/32.0;
prolongation[1](4,3) = -35.0/64.0;
prolongation[1](4,4) = 35.0/32.0;
+
+ restriction[0](0,0) = 1.0;
+ restriction[0](2,3) = 1.0;
+ restriction[0](3,1) = 1.0;
+ restriction[1](1,1) = 1.0;
+ restriction[1](3,0) = 1.0;
+ restriction[1](4,3) = 1.0;
};
prolongation[3](23,19) = 35.0/32.0;
prolongation[3](23,22) = -35.0/64.0;
prolongation[3](24,19) = 1.0;
+
+ restriction[0](0,0) = 1.0;
+ restriction[0](4,5) = 1.0;
+ restriction[0](5,1) = 1.0;
+ restriction[0](13,14) = 1.0;
+ restriction[0](14,3) = 1.0;
+ restriction[0](16,24) = 1.0;
+ restriction[0](20,8) = 1.0;
+ restriction[0](23,11) = 1.0;
+ restriction[0](24,2) = 1.0;
+ restriction[1](1,1) = 1.0;
+ restriction[1](5,0) = 1.0;
+ restriction[1](6,5) = 1.0;
+ restriction[1](7,8) = 1.0;
+ restriction[1](8,2) = 1.0;
+ restriction[1](17,24) = 1.0;
+ restriction[1](20,14) = 1.0;
+ restriction[1](21,11) = 1.0;
+ restriction[1](24,3) = 1.0;
+ restriction[2](2,2) = 1.0;
+ restriction[2](8,1) = 1.0;
+ restriction[2](9,8) = 1.0;
+ restriction[2](11,3) = 1.0;
+ restriction[2](12,11) = 1.0;
+ restriction[2](18,24) = 1.0;
+ restriction[2](21,5) = 1.0;
+ restriction[2](22,14) = 1.0;
+ restriction[2](24,0) = 1.0;
+ restriction[3](3,3) = 1.0;
+ restriction[3](10,11) = 1.0;
+ restriction[3](11,2) = 1.0;
+ restriction[3](14,0) = 1.0;
+ restriction[3](15,14) = 1.0;
+ restriction[3](19,24) = 1.0;
+ restriction[3](22,8) = 1.0;
+ restriction[3](23,5) = 1.0;
+ restriction[3](24,1) = 1.0;
};