From bcbfae7861c069c4775363a561738a2766e77284 Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 3 Nov 1998 15:53:21 +0000 Subject: [PATCH] Implement the restriction (interpolation) matrices for transfer from child to parent cells. git-svn-id: https://svn.dealii.org/trunk@622 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe.h | 26 +- .../deal.II/source/fe/fe_lib.criss_cross.cc | 54 ++++ deal.II/deal.II/source/fe/fe_lib.cubic.cc | 127 ++++++++ deal.II/deal.II/source/fe/fe_lib.linear.cc | 40 --- deal.II/deal.II/source/fe/fe_lib.quadratic.cc | 294 +++++++++++++++++- deal.II/deal.II/source/fe/fe_lib.quartic.cc | 121 +++++++ 6 files changed, 620 insertions(+), 42 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 29d3f03701..f66146930a 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -236,7 +236,8 @@ struct FiniteElementBase : * 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; @@ -309,6 +310,29 @@ struct FiniteElementBase : * 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 diff --git a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc index 793a902255..8d984409dd 100644 --- a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc +++ b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc @@ -62,6 +62,50 @@ 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. @@ -238,6 +282,7 @@ 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): @@ -442,6 +487,15 @@ FECrissCross<2>::FECrissCross () : 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; }; diff --git a/deal.II/deal.II/source/fe/fe_lib.cubic.cc b/deal.II/deal.II/source/fe/fe_lib.cubic.cc index 5529d90e5a..60018e382c 100644 --- a/deal.II/deal.II/source/fe/fe_lib.cubic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.cubic.cc @@ -68,6 +68,34 @@ 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, @@ -79,6 +107,7 @@ 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); ----------------------------------------------------------------------- @@ -89,6 +118,7 @@ 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 */ @@ -211,6 +241,50 @@ 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 @@ -277,6 +351,7 @@ 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); @@ -293,6 +368,8 @@ 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 */ @@ -338,6 +415,39 @@ FECubicSub<1>::FECubicSub () : 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; }; @@ -867,6 +977,23 @@ FECubicSub<2>::FECubicSub () : 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; }; diff --git a/deal.II/deal.II/source/fe/fe_lib.linear.cc b/deal.II/deal.II/source/fe/fe_lib.linear.cc index 4336f78321..e8c18dcf3e 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -31,11 +31,6 @@ FELinear<1>::FELinear () : // 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; @@ -138,44 +133,9 @@ FELinear<2>::FELinear () : 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.; diff --git a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc index 793f36d5b8..e0fdf45f74 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc @@ -11,6 +11,244 @@ +/*----------------------------------------------------------------- + * 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 @@ -29,6 +267,7 @@ FEQuadraticSub<1>::FEQuadraticSub () : 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 @@ -38,8 +277,39 @@ FEQuadraticSub<1>::FEQuadraticSub () : 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; @@ -60,6 +330,11 @@ FEQuadraticSub<1>::FEQuadraticSub () : 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; }; @@ -163,7 +438,7 @@ FEQuadraticSub<2>::FEQuadraticSub () : 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; @@ -525,6 +800,23 @@ FEQuadraticSub<2>::FEQuadraticSub () : 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; }; diff --git a/deal.II/deal.II/source/fe/fe_lib.quartic.cc b/deal.II/deal.II/source/fe/fe_lib.quartic.cc index 50ba7da68e..8170590e5c 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quartic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quartic.cc @@ -76,10 +76,38 @@ 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); ----------------------------------------------------------------------- @@ -93,6 +121,8 @@ 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 */ @@ -228,6 +258,50 @@ 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 @@ -293,6 +367,7 @@ 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); @@ -310,6 +385,8 @@ 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 */ @@ -349,6 +426,13 @@ FEQuarticSub<1>::FEQuarticSub () : 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; }; @@ -1176,6 +1260,43 @@ FEQuarticSub<2>::FEQuarticSub () : 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; }; -- 2.39.5