]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement the restriction (interpolation) matrices for transfer from child to parent...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 3 Nov 1998 15:53:21 +0000 (15:53 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 3 Nov 1998 15:53:21 +0000 (15:53 +0000)
git-svn-id: https://svn.dealii.org/trunk@622 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/source/fe/fe_lib.criss_cross.cc
deal.II/deal.II/source/fe/fe_lib.cubic.cc
deal.II/deal.II/source/fe/fe_lib.linear.cc
deal.II/deal.II/source/fe/fe_lib.quadratic.cc
deal.II/deal.II/source/fe/fe_lib.quartic.cc

index 29d3f03701d60959abf7c8f1262c5ebb64eaa78c..f66146930a2273750b21ee3b90f04ed0ac8c1524 100644 (file)
@@ -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
index 793a9022552721a446776a0aee549cdd28fe6823..8d984409dd2555db529622395562335b389f72ac 100644 (file)
     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):
   
@@ -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;
 };
 
 
index 5529d90e5ab920844b4ead979b6ec13bd060455b..60018e382c18d2a1aa52bb84ed049b46bd5e86c3 100644 (file)
     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
 */
 
@@ -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;
 };
 
 
index 4336f783215620cc0920479cfd75301ff03e25b1..e8c18dcf3e083153ee4fc5db1ddfa6cb65b1e925 100644 (file)
@@ -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.;
index 793f36d5b8f9bd7b90516145b08b8a3a774f2880..e0fdf45f743399313ae533f0c6b4c3263aca369e 100644 (file)
 
 
 
+/*-----------------------------------------------------------------
+ * 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; 
 };
 
 
index 50ba7da68e64a36e54018fdc6041759f3d87d19d..8170590e5c0ce211196db7eedff4374c5833a696 100644 (file)
     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
 */
 
@@ -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;
 };
 
 

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.