* 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}
*
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()));
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
-/*--------------------------------- 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
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);
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);
-/*-----------------------------------------------------------------
- * 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 <>
-/*--------------------------------- 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
-*/
-
-
-
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());
};
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());
};
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());
};
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()));
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);
};
};
+
+ 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);
};
(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)