]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Invent a faster scheme to compute the jacobian matrix from the vertices and the point...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 4 Nov 1998 18:31:44 +0000 (18:31 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 4 Nov 1998 18:31:44 +0000 (18:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@631 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_values.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.quadratic.cc
deal.II/deal.II/source/fe/fe_lib.quartic.cc
deal.II/deal.II/source/fe/fe_linear_mapping.cc
deal.II/deal.II/source/fe/fe_values.cc

index 6f319388d3a3aec819551c0261cb4dc5fb265eb9..35ce803a3a0065e5cbeb7f214b26e3f0f2492054 100644 (file)
@@ -107,6 +107,17 @@ template <int dim> class Quadrature;
  *  transformation from unit to real face to compute the determinant of the
  *  Jacobi matrix to get the scaling of the surface element $do$.
  *
+ *  The question whether to compute the Jacobi matrix as the inverse of another
+ *  matrix M (which we can compute from the transformation, while we can't do
+ *  so for the Jacobi matrix itself) or its transpose is a bit delicate. It
+ *  should be kept in mind that when we compute the gradients in real space
+ *  from those on the unit cell, we multiply with the Jacobi matrix
+ *  \textit{from the right}; the whole situation is a bit confusing and it
+ *  either takes deep though or trial-and-error to do it right. Some more
+ *  information on this can be found in the source code documentation for the
+ *  #FELinearMapping<dim>::fill_fe_values# function, where also a small test
+ *  program is presented.
+ *
  *  
  *  \subsection{Member functions}
  *
index 8d984409dd2555db529622395562335b389f72ac..7244c61d43dc1fb8c53e7c4bd5ffd1065e5fe584 100644 (file)
@@ -911,7 +911,7 @@ void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &ce
                                        vector<Point<dim> > &q_points,
                                        const bool           compute_q_points,
                                        const dFMatrix      &shape_values_transform,
-                                       const vector<vector<Point<dim> > > &shape_grad_transform,
+                                       const vector<vector<Point<dim> > > &/*shape_grad_transform*/,
                                        const Boundary<dim> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
@@ -964,25 +964,53 @@ void FECrissCross<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator &ce
 
   However, we rewrite the loops to only compute the gradient once for
   each integration point and basis function.
+
+  The scheme laid down above was originally used. Due to recent advances
+  in the authors understanding of most basic things, it was dropped and
+  replaced by the following version. See #FELinearMapping<dim>::fill_fe_values#
+  for more information on this.
 */
-  if (compute_jacobians) 
-    {
-      dFMatrix M(dim,dim);
-      for (unsigned int l=0; l<n_points; ++l) 
-       {
-         M.clear ();
-         for (unsigned int s=0; s<n_transform_functions; ++s)
-           {
-                                              // we want the linear transform,
-                                              // so use that function
-             const Point<dim> gradient = shape_grad_transform[s][l];
-             for (unsigned int i=0; i<dim; ++i)
-               for (unsigned int j=0; j<dim; ++j)
-                 M(i,j) += support_points[s](i) * gradient(j);
-           };
-         jacobians[l].invert(M);
-       };
-    };
+
+  Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
+  for (unsigned int l=0; l<GeometryInfo<dim>::vertices_per_cell; ++l)
+    vertices[l] = cell->vertex(l);
+
+  if (compute_jacobians)
+    for (unsigned int point=0; point<n_points; ++point)
+      {
+       const double xi = unit_points[point](0);
+       const double eta= unit_points[point](1);
+       
+       const double t6 = vertices[0](0)*vertices[3](1);
+       const double t8 = vertices[2](0)*xi;
+       const double t10 = vertices[1](0)*eta;
+       const double t12 = vertices[3](0)*vertices[1](1);
+       const double t16 = vertices[3](0)*xi;
+       const double t20 = vertices[0](0)*vertices[1](1);
+       const double t22 = vertices[0](0)*vertices[2](1);
+       const double t24 = t6*xi-t8*vertices[1](1)-t10*vertices[3](1)+
+                          t12*eta-vertices[3](0)*vertices[2](1)*eta-
+                          t16*vertices[0](1)+t16*vertices[1](1)-t12+
+                          vertices[3](0)*vertices[0](1)-t20*eta+t22*eta;
+       const double t28 = vertices[1](0)*vertices[3](1);
+       const double t31 = vertices[2](0)*eta;
+       const double t36 = t8*vertices[0](1)+vertices[1](0)*vertices[2](1)*xi-
+                          t28*xi+t10*vertices[0](1)-t31*vertices[0](1)+
+                          t31*vertices[3](1)+t20-t6-vertices[1](0)*
+                          vertices[0](1)+t28-t22*xi;
+       const double t38 = 1/(t24+t36);
+
+       jacobians[point](0,0) = (-vertices[0](1)+vertices[0](1)*xi-
+                                vertices[1](1)*xi+vertices[2](1)*xi+
+                                vertices[3](1)-vertices[3](1)*xi)*t38;
+       jacobians[point](0,1) = -(-vertices[0](0)+vertices[0](0)*xi-
+                                 vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38;
+       jacobians[point](1,0) = -(-vertices[0](1)+vertices[0](1)*eta+
+                                 vertices[1](1)-vertices[1](1)*eta+
+                                 vertices[2](1)*eta-vertices[3](1)*eta)*t38;
+       jacobians[point](1,1) = (-vertices[0](0)+vertices[0](0)*eta+
+                                vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38;
+      };
 };
 
 #endif
index 60018e382c18d2a1aa52bb84ed049b46bd5e86c3..71c5b35670dd3a3254145b1bcccc8d3bae5f704d 100644 (file)
 
 
 
-/*--------------------------------- 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.
-  --------------------------------------------------------------------------
-  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:
-  
-  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
@@ -454,8 +203,8 @@ FECubicSub<1>::FECubicSub () :
 
 template <>
 double
-FECubicSub<1>::shape_value(const unsigned int i,
-                              const Point<1>     &p) const
+FECubicSub<1>::shape_value (const unsigned int i,
+                           const Point<1>     &p) const
 {
   Assert((i<total_dofs), ExcInvalidIndex(i));
   const double xi = p(0);
@@ -473,8 +222,8 @@ FECubicSub<1>::shape_value(const unsigned int i,
 
 template <>
 Point<1>
-FECubicSub<1>::shape_grad(const unsigned int i,
-                             const Point<1>    &p) const
+FECubicSub<1>::shape_grad (const unsigned int i,
+                          const Point<1>    &p) const
 {
   Assert((i<total_dofs), ExcInvalidIndex(i));
   const double xi = p(0);
index e0fdf45f743399313ae533f0c6b4c3263aca369e..5605d05b2873ce24c5a0ab807c5f3ad13aca380e 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
 
 template <>
index 8170590e5c0ce211196db7eedff4374c5833a696..1bcda503d5bd2122faa95dcf1893b6a491e4e27d 100644 (file)
 
 
 
-/*--------------------------------- 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.
-  --------------------------------------------------------------------------
-  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:
-  
-  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
-*/
-
-
-
 
 
 
index 5e1b2cb8b8545c712485275f6fa00164a6445709..bc1fda90ea640ce74158ccc0238d2008dc27fff8 100644 (file)
@@ -53,9 +53,9 @@ FELinearMapping<1>::shape_grad_transform(const unsigned int i,
 
 template <>
 void FELinearMapping<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &,
-                                           const Boundary<1>         &,
-                                           const vector<Point<0> > &,
-                                           vector<double>      &) const {
+                                            const Boundary<1>         &,
+                                            const vector<Point<0> > &,
+                                            vector<double>      &) const {
   Assert (false, ExcInternalError());
 };
 
@@ -63,9 +63,9 @@ void FELinearMapping<1>::get_face_jacobians (const DoFHandler<1>::face_iterator
 
 template <>
 void FELinearMapping<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &,
-                                              const unsigned int           ,
-                                              const vector<Point<0> > &,
-                                              vector<double>      &) const {
+                                               const unsigned int           ,
+                                               const vector<Point<0> > &,
+                                               vector<double>      &) const {
   Assert (false, ExcInternalError());
 };
 
@@ -73,10 +73,10 @@ void FELinearMapping<1>::get_subface_jacobians (const DoFHandler<1>::face_iterat
 
 template <>
 void FELinearMapping<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &,
-                                           const unsigned int,
-                                           const Boundary<1> &,
-                                           const vector<Point<0> > &,
-                                           vector<Point<1> > &) const {
+                                            const unsigned int,
+                                            const Boundary<1> &,
+                                            const vector<Point<0> > &,
+                                            vector<Point<1> > &) const {
   Assert (false, ExcInternalError());
 };
 
@@ -286,7 +286,7 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
                                           vector<Point<dim> > &q_points,
                                           const bool           compute_q_points,
                                           const dFMatrix      &shape_values_transform,
-                                          const vector<vector<Point<dim> > > &shape_grad_transform,
+                                          const vector<vector<Point<dim> > > &/*shape_grad_transform*/,
                                           const Boundary<dim> &boundary) const {
   Assert (jacobians.size() == unit_points.size(),
          ExcWrongFieldDimension(jacobians.size(), unit_points.size()));
@@ -344,7 +344,13 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
 
   However, we rewrite the loops to only compute the gradient once for
   each integration point and basis function.
-*/
+
+  Indeed, this was the old way we did it; the code is below. However, there
+  is a more efficient way, namely setting up M analytically, doing the
+  inversion analyically and only then doing the evaluation; a small Maple
+  script does this (it is part of the <lagrange> script in the <scripts>
+  subdirectory).
+  
   if (compute_jacobians) 
     {
       dFMatrix M(dim,dim);
@@ -364,6 +370,101 @@ void FELinearMapping<dim>::fill_fe_values (const DoFHandler<dim>::cell_iterator
        };
     };
 
+    
+  One last note regarding whether we have to invert M or M transposed: it is
+  easy to try out, by computing the gradients of a function on a distorted
+  cell (just move one vertex) where the nodal values for linear elements
+  are one for the moved vertex and zero otherwise. Please also note that
+  when computing the gradients on the real cell, the jacobian matrix
+  is multiplied to the unit cell gradient *from the right*! be very careful
+  with these things.
+
+  The following little program tests the correct behaviour (you have to find
+  out the right include files, I tested it within a whole project with far
+  more include files than necessary):
+
+  -------------------------------------------
+  int main () {
+    Triangulation<2> tria;
+    tria.create_hypercube (0,1);
+    tria.begin_active()->vertex(2)(0) = 2;
+
+    DoFHandler<2> dof(&tria);
+    FELinear<2> fe;
+    dof.distribute_dofs(fe);
+
+    StraightBoundary<2> b;
+    QTrapez<2> q;
+    FEValues<2> fevalues(fe,q,update_gradients);
+    fevalues.reinit (dof.begin_active(),b);
+  
+  
+    dVector val(4);
+    val(2) = 1;
+
+    vector<Point<2> > grads(4);
+    fevalues.get_function_grads (val, grads);
+
+    for (unsigned int i=0; i<4; ++i)
+      cout << "Vertex " << i
+          << "   grad=" << grads[i] << endl;
+  };
+  ---------------------------------------------
+  
+  The correct output should be
+  --------------------------------
+  Vertex 0   grad=0 0
+  Vertex 1   grad=0.5 0
+  Vertex 2   grad=0 1
+  Vertex 3   grad=0.5 0.5
+  --------------------------------
+  and the wrong would be
+  --------------------------------
+  Vertex 0   grad=0 0
+  Vertex 1   grad=0.5 0
+  Vertex 2   grad=-1 1
+  Vertex 3   grad=0 1
+  --------------------------------  
+*/
+
+  if (compute_jacobians)
+    for (unsigned int point=0; point<n_points; ++point)
+      {
+       const double xi = unit_points[point](0);
+       const double eta= unit_points[point](1);
+       
+       const double t6 = vertices[0](0)*vertices[3](1);
+       const double t8 = vertices[2](0)*xi;
+       const double t10 = vertices[1](0)*eta;
+       const double t12 = vertices[3](0)*vertices[1](1);
+       const double t16 = vertices[3](0)*xi;
+       const double t20 = vertices[0](0)*vertices[1](1);
+       const double t22 = vertices[0](0)*vertices[2](1);
+       const double t24 = t6*xi-t8*vertices[1](1)-t10*vertices[3](1)+
+                          t12*eta-vertices[3](0)*vertices[2](1)*eta-
+                          t16*vertices[0](1)+t16*vertices[1](1)-t12+
+                          vertices[3](0)*vertices[0](1)-t20*eta+t22*eta;
+       const double t28 = vertices[1](0)*vertices[3](1);
+       const double t31 = vertices[2](0)*eta;
+       const double t36 = t8*vertices[0](1)+vertices[1](0)*vertices[2](1)*xi-
+                          t28*xi+t10*vertices[0](1)-t31*vertices[0](1)+
+                          t31*vertices[3](1)+t20-t6-vertices[1](0)*
+                          vertices[0](1)+t28-t22*xi;
+       const double t38 = 1/(t24+t36);
+
+       jacobians[point](0,0) = (-vertices[0](1)+vertices[0](1)*xi-
+                                vertices[1](1)*xi+vertices[2](1)*xi+
+                                vertices[3](1)-vertices[3](1)*xi)*t38;
+       jacobians[point](0,1) = -(-vertices[0](0)+vertices[0](0)*xi-
+                                 vertices[1](0)*xi+t8+vertices[3](0)-t16)*t38;
+       jacobians[point](1,0) = -(-vertices[0](1)+vertices[0](1)*eta+
+                                 vertices[1](1)-vertices[1](1)*eta+
+                                 vertices[2](1)*eta-vertices[3](1)*eta)*t38;
+       jacobians[point](1,1) = (-vertices[0](0)+vertices[0](0)*eta+
+                                vertices[1](0)-t10+t31-vertices[3](0)*eta)*t38;
+      };
+  
+    
   if (compute_support_points)
     get_support_points (cell, boundary, support_points);
 };
index 89a9141faedd24ad364070085a805cefd48c1062..3957483e7a1056c1d8ebaf4391bb6f95b1b6b270 100644 (file)
@@ -214,18 +214,18 @@ void FEValues<dim>::reinit (const typename DoFHandler<dim>::cell_iterator &cell,
       (update_flags & update_gradients) ||
       (update_flags & update_support_points))
     fe->fill_fe_values (cell,
-                      unit_quadrature_points,
-                      jacobi_matrices,
-                      update_flags & (update_jacobians  |
-                                      update_JxW_values |
-                                      update_gradients),
-                      support_points,
-                      update_flags & update_support_points,
-                      quadrature_points,
-                      update_flags & update_q_points,
-                      shape_values_transform[0], unit_shape_gradients_transform,
-                      boundary);
-
+                       unit_quadrature_points,
+                       jacobi_matrices,
+                       update_flags & (update_jacobians  |
+                                       update_JxW_values |
+                                       update_gradients),
+                       support_points,
+                       update_flags & update_support_points,
+                       quadrature_points,
+                       update_flags & update_q_points,
+                       shape_values_transform[0], unit_shape_gradients_transform,
+                       boundary);
+  
                                   // compute gradients on real element if
                                   // requested
   if (update_flags & update_gradients) 

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.