From d0778b40302ff631c6dd09f2fbd00cfe2ff6979c Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 4 Nov 1998 18:31:44 +0000 Subject: [PATCH] Invent a faster scheme to compute the jacobian matrix from the vertices and the point of integration. Add some doc. git-svn-id: https://svn.dealii.org/trunk@631 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_values.h | 11 + .../deal.II/source/fe/fe_lib.criss_cross.cc | 66 +++-- deal.II/deal.II/source/fe/fe_lib.cubic.cc | 259 +---------------- deal.II/deal.II/source/fe/fe_lib.quadratic.cc | 240 ---------------- deal.II/deal.II/source/fe/fe_lib.quartic.cc | 263 ------------------ .../deal.II/source/fe/fe_linear_mapping.cc | 125 ++++++++- deal.II/deal.II/source/fe/fe_values.cc | 24 +- 7 files changed, 187 insertions(+), 801 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index 6f319388d3..35ce803a3a 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -107,6 +107,17 @@ template 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::fill_fe_values# function, where also a small test + * program is presented. + * * * \subsection{Member functions} * diff --git a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc index 8d984409dd..7244c61d43 100644 --- a/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc +++ b/deal.II/deal.II/source/fe/fe_lib.criss_cross.cc @@ -911,7 +911,7 @@ void FECrissCross::fill_fe_values (const DoFHandler::cell_iterator &ce vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, - const vector > > &shape_grad_transform, + const vector > > &/*shape_grad_transform*/, const Boundary &boundary) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); @@ -964,25 +964,53 @@ void FECrissCross::fill_fe_values (const DoFHandler::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::fill_fe_values# + for more information on this. */ - if (compute_jacobians) - { - dFMatrix M(dim,dim); - for (unsigned int l=0; l gradient = shape_grad_transform[s][l]; - for (unsigned int i=0; i vertices[GeometryInfo::vertices_per_cell]; + for (unsigned int l=0; l::vertices_per_cell; ++l) + vertices[l] = cell->vertex(l); + + if (compute_jacobians) + for (unsigned int point=0; point1/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::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 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 <> diff --git a/deal.II/deal.II/source/fe/fe_lib.quartic.cc b/deal.II/deal.II/source/fe/fe_lib.quartic.cc index 8170590e5c..1bcda503d5 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quartic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quartic.cc @@ -129,269 +129,6 @@ -/*--------------------------------- 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 -*/ - - - diff --git a/deal.II/deal.II/source/fe/fe_linear_mapping.cc b/deal.II/deal.II/source/fe/fe_linear_mapping.cc index 5e1b2cb8b8..bc1fda90ea 100644 --- a/deal.II/deal.II/source/fe/fe_linear_mapping.cc +++ b/deal.II/deal.II/source/fe/fe_linear_mapping.cc @@ -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 > &, - vector &) const { + const Boundary<1> &, + const vector > &, + vector &) 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 > &, - vector &) const { + const unsigned int , + const vector > &, + vector &) 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 > &, - vector > &) const { + const unsigned int, + const Boundary<1> &, + const vector > &, + vector > &) const { Assert (false, ExcInternalError()); }; @@ -286,7 +286,7 @@ void FELinearMapping::fill_fe_values (const DoFHandler::cell_iterator vector > &q_points, const bool compute_q_points, const dFMatrix &shape_values_transform, - const vector > > &shape_grad_transform, + const vector > > &/*shape_grad_transform*/, const Boundary &boundary) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); @@ -344,7 +344,13 @@ void FELinearMapping::fill_fe_values (const DoFHandler::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 script in the + subdirectory). + if (compute_jacobians) { dFMatrix M(dim,dim); @@ -364,6 +370,101 @@ void FELinearMapping::fill_fe_values (const DoFHandler::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 > 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::reinit (const typename DoFHandler::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) -- 2.39.5