From 9012ddd7b25ed68a8f3804a0113ac8abf12f111a Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 25 Jun 1998 15:00:42 +0000 Subject: [PATCH] Some further changes for (bi-)quadratic elements. git-svn-id: https://svn.dealii.org/trunk@412 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_lib.linear.cc | 362 ++++++++++++++++++--- 1 file changed, 324 insertions(+), 38 deletions(-) diff --git a/deal.II/deal.II/source/fe/fe_lib.linear.cc b/deal.II/deal.II/source/fe/fe_lib.linear.cc index cfdf0bb204..11ef7c4980 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -766,6 +766,28 @@ void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator & Assert (false, ExcInternalError()); }; + + +template <> +void FEQuadraticSub<1>::get_local_mass_matrix (const DoFHandler<1>::cell_iterator &cell, + const Boundary<1> &, + dFMatrix &local_mass_matrix) const { + Assert (local_mass_matrix.n() == total_dofs, + ExcWrongFieldDimension(local_mass_matrix.n(),total_dofs)); + Assert (local_mass_matrix.m() == total_dofs, + ExcWrongFieldDimension(local_mass_matrix.m(),total_dofs)); + + const double h = cell->vertex(1)(0) - cell->vertex(0)(0); + Assert (h>0, ExcJacobiDeterminantHasWrongSign()); + + local_mass_matrix(0,0) = local_mass_matrix(1,1) = 2./15.*h; + local_mass_matrix(0,1) = local_mass_matrix(1,0) = -1./30.*h; + local_mass_matrix(0,2) = local_mass_matrix(2,0) = 1./15.*h; + local_mass_matrix(1,2) = local_mass_matrix(2,1) = 1./15.*h; + local_mass_matrix(2,2) = local_mass_matrix(2,2) = 8./15.*h; +}; + + #endif @@ -800,14 +822,14 @@ FEQuadraticSub<2>::shape_value (const unsigned int i, eta= p(1); switch (i) { - case 0: return (1.-xi)*( 2*xi-1) * (1.-eta)*( 2*eta-1); - case 1: return xi *(-2*xi+1) * (1.-eta)*( 2*eta-1); - case 2: return xi *(-2*xi+1) * eta *(-2*eta+1); - case 3: return (1.-xi)*( 2*xi-1) * eta *(-2*eta+1); - case 4: return 4 * (1-xi)*xi * (1-eta)*(1-2*eta); + case 0: return (1-xi)*( 2*xi-1) * (1-eta)*( 2*eta-1); + case 1: return xi *(-2*xi+1) * (1-eta)*( 2*eta-1); + case 2: return xi *(-2*xi+1) * eta *(-2*eta+1); + case 3: return (1-xi)*( 2*xi-1) * eta *(-2*eta+1); + case 4: return 4 * (1-xi)*xi * (1-eta)*(1-2*eta); case 5: return 4 * xi *(-1+2*xi) * (1-eta)*eta; - case 6: return 4 * (1-xi)*xi * eta *(-1+2*eta); - case 7: return 4 * (1.-xi)*(1-2*xi) * (1-eta)*eta; + case 6: return 4 * (1-xi)*xi * eta *(-1+2*eta); + case 7: return 4 * (1-xi)*(1-2*xi) * (1-eta)*eta; case 8: return 16 * xi*(1-xi) * eta*(1-eta); }; return 0; @@ -826,29 +848,306 @@ FEQuadraticSub<2>::shape_grad (const unsigned int i, eta= p(1); switch (i) { - case 0: return Point<2>((-4*xi+3) * (1.-eta)*( 2*eta-1), - (1.-xi)*( 2*xi-1) * (-4*eta+3)); - case 1: return Point<2>((-4*xi+1) * (1.-eta)*( 2*eta-1) , - xi *(-2*xi+1) * (-4*eta+3)); - case 2: return Point<2>((-4*xi+1) * eta *(-2*eta+1), - xi *(-2*xi+1) * (-4*eta+1)); - case 3: return Point<2>((-4*xi+3) * eta *(-2*eta+1), - (1.-xi)*( 2*xi-1) * (-4*eta+1)); - case 4: return Point<2>(4 * (1-2*xi) * (1-eta)*(1-2*eta), - 4 * (1-xi)*xi * (4*eta-3)); - case 5: return Point<2>(4 * (4*xi-1) * (1-eta)*eta, - 4 * xi *(-1+2*xi) * (1-2*eta)); - case 6: return Point<2>(4 * (1-2*xi) * eta *(-1+2*eta), - 4 * (1-xi)*xi * (4*eta-1)); - case 7: return Point<2>(4 * (4*xi-3) * (1-eta)*eta, - 4 * (1.-xi)*(1-2*xi) * (1-2*eta)); - case 8: return Point<2>(16 * (1-2*xi) * eta*(1-eta), - 16 * xi*(1-xi) * (1-2*eta)); + case 0: return Point<2>(-(2*xi-1)*(1-eta)*(2*eta-1)+2*(1-xi)*(1-eta)*(2*eta-1), + -(1-xi)*(2*xi-1)*(2*eta-1)+2*(1-xi)*(2*xi-1)*(1-eta)); + case 1: return Point<2>((-2*xi+1)*(1-eta)*(2*eta-1)-2*xi*(1-eta)*(2*eta-1), + -xi*(-2*xi+1)*(2*eta-1)+2*xi*(-2*xi+1)*(1-eta)); + case 2: return Point<2>((-2*xi+1)*eta*(-2*eta+1)-2*xi*eta*(-2*eta+1), + xi*(-2*xi+1)*(-2*eta+1)-2*xi*(-2*xi+1)*eta); + case 3: return Point<2>(-(2*xi-1)*eta*(-2*eta+1)+2*(1-xi)*eta*(-2*eta+1), + (1-xi)*(2*xi-1)*(-2*eta+1)-2*(1-xi)*(2*xi-1)*eta); + case 4: return Point<2>(-4*xi*(1-eta)*(-2*eta+1)+4*(1-xi)*(1-eta)*(-2*eta+1), + -4*(1-xi)*xi*(-2*eta+1)-8*(1-xi)*xi*(1-eta)); + case 5: return Point<2>(4*(2*xi-1)*(1-eta)*eta+8*xi*(1-eta)*eta, + -4*xi*(2*xi-1)*eta+4*xi*(2*xi-1)*(1-eta)); + case 6: return Point<2>(-4*xi*eta*(2*eta-1)+4*(1-xi)*eta*(2*eta-1), + 4*(1-xi)*xi*(2*eta-1)+8*(1-xi)*xi*eta); + case 7: return Point<2>(-4*(-2*xi+1)*(1-eta)*eta-8*(1-xi)*(1-eta)*eta, + -4*(1-xi)*(-2*xi+1)*eta+4*(1-xi)*(-2*xi+1)*(1-eta)); + case 8: return Point<2>(16*(1-xi)*(1-eta)*eta-16*xi*eta*(1-eta), + 16*xi*(1-xi)*(1-eta)-16*(1-xi)*xi*eta); }; return Point<2> (); }; + +template <> +void FEQuadraticSub<2>::get_local_mass_matrix (const DoFHandler<2>::cell_iterator &cell, + const Boundary<2> &, + dFMatrix &local_mass_matrix) const { + Assert (local_mass_matrix.n() == total_dofs, + ExcWrongFieldDimension(local_mass_matrix.n(),total_dofs)); + Assert (local_mass_matrix.m() == total_dofs, + ExcWrongFieldDimension(local_mass_matrix.m(),total_dofs)); + +/* Get the computation of the local mass matrix by these lines in maple. Note + that tphi[i] are the basis function of the linear finite element, which + are used by the transformation (therefore >tvertex(0)(0), + cell->vertex(1)(0), + cell->vertex(2)(0), + cell->vertex(3)(0) }; + const double y[4] = { cell->vertex(0)(1), + cell->vertex(1)(1), + cell->vertex(2)(1), + cell->vertex(3)(1) }; + +/* check that the Jacobi determinant + + t0 = (-x[0]*(1.0-eta)+x[1]*(1.0-eta)+x[2]*eta-x[3]*eta) * + (-y[0]*(1.0-xi)-y[1]*xi+y[2]*xi+y[3]*(1.0-xi)) - + (-x[0]*(1.0-xi)-x[1]*xi+x[2]*xi+x[3]*(1.0-xi)) * + (-y[0]*(1.0-eta)+y[1]*(1.0-eta)+y[2]*eta-y[3]*eta) + + has the right sign. + + We do not attempt to check its (hopefully) positive sign at all points + on the unit cell, but we check that it is positive in the four corners, + which is sufficient since $det J$ is a bilinear function. +*/ + Assert ((-x[0]+x[1])*(-y[0]+y[3])-(-x[0]+x[3])*(-y[0]+y[1]), // xi=eta=0 + ExcJacobiDeterminantHasWrongSign()); + Assert ((x[2]-x[3])*(-y[0]+y[3])-(-x[0]+x[3])*(y[2]-y[3]), // xi=0, eta=1 + ExcJacobiDeterminantHasWrongSign()); + Assert ((x[2]-x[3])*(-y[1]+y[2])-(-x[1]+x[2])*(y[2]-y[3]), // xi=eta=1 + ExcJacobiDeterminantHasWrongSign()); + Assert ((-x[0]+x[1])*(-y[1]+y[2])-(-x[1]+x[2])*(-y[0]+y[1]), // xi=1, eta=0 + ExcJacobiDeterminantHasWrongSign()); + + const double t1 = (x[1]*y[0]); + const double t2 = (x[1]*y[2]); + const double t3 = (x[0]*y[3]); + const double t4 = (x[3]*y[2]); + const double t5 = (x[2]*y[3]); + const double t6 = (x[0]*y[1]); + const double t7 = (x[3]*y[1]); + const double t8 = (x[3]*y[0]); + const double t9 = (x[2]*y[1]); + const double t10 = (x[1]*y[3]); + const double t12 = (x[0]*y[2]); + const double t13 = (x[2]*y[0]); + const double t14 = (7.0/1800.0*t1-t2/450+t3/450+t4/1800-t5/1800- + 7.0/1800.0*t6+t12/600+ + t7/600-t8/450-t13/600+t9/450-t10/600); + const double t15 = (-t1/1800+t2/1800-t3/1800-t4/1800+t5/1800+ + t6/1800+t8/1800-t9/1800); + const double t16 = (t1/450-t2/1800+7.0/1800.0*t3+t4/450- + t5/450-t6/450-t12/600+t7/600 + -7.0/1800.0*t8+t13/600+t9/1800-t10/600); + const double t17 = (-7.0/900.0*t1-2.0/225.0*t3-t4/900+t5/900 + +7.0/900.0*t6+t12/900-7.0/ + 900.0*t7+2.0/225.0*t8-t13/900+7.0/900.0*t10); + const double t18 = (t1/450-t2/900+t3/900-t6/450+t12/900+ + t7/900-t8/900-t13/900+t9/900- + t10/900); + const double t19 = (t1/900+t3/450+t4/900-t5/900-t6/900 + -t12/900+t7/900-t8/450+t13/900- + t10/900); + const double t20 = (-2.0/225.0*t1+t2/900-7.0/900.0*t3+ + 2.0/225.0*t6-t12/900-7.0/900.0*t7 + +7.0/900.0*t8+t13/900-t9/900+7.0/900.0*t10); + const double t21 = (-t1/225-t3/225+t6/225-t7/225+t8/225+t10/225); + const double t23 = (t1/450-7.0/1800.0*t2+t3/1800+t4/450 + -t5/450-t6/450+t12/600-t7/600-t8 + /1800-t13/600+7.0/1800.0*t9+t10/600); + const double t24 = (-7.0/900.0*t1+2.0/225.0*t2-t4/900+t5/900 + +7.0/900.0*t6-7.0/900.0*t12 + +t7/900+7.0/900.0*t13-2.0/225.0*t9-t10/900); + const double t25 = (-2.0/225.0*t1+7.0/900.0*t2-t3/900+2.0/225.0*t6 + -7.0/900.0*t12-t7/900 + +t8/900+7.0/900.0*t13-7.0/900.0*t9+t10/900); + const double t26 = (t1/900-t2/450+t4/900-t5/900-t6/900+t12/900 + -t7/900-t13/900+t9/450+ + t10/900); + const double t27 = (-t1/225+t2/225+t6/225-t12/225+t13/225-t9/225); + const double t29 = (t1/1800-t2/450+t3/450+7.0/1800.0*t4-7.0/1800.0*t5 + -t6/1800-t12/600- + t7/600-t8/450+t13/600+t9/450+t10/600); + const double t30 = (7.0/900.0*t2-t3/900-2.0/225.0*t4+2.0/225.0*t5 + +t12/900+7.0/900.0*t7+ + t8/900-t13/900-7.0/900.0*t9-7.0/900.0*t10); + const double t31 = (-t1/900+2.0/225.0*t2-7.0/900.0*t4+7.0/900.0*t5 + +t6/900-t12/900+7.0/ + 900.0*t7+t13/900-2.0/225.0*t9-7.0/900.0*t10); + const double t32 = (-t2/900+t3/900+t4/450-t5/450-t12/900-t7/900 + -t8/900+t13/900+t9/900+ + t10/900); + const double t33 = (t2/225-t4/225+t5/225+t7/225-t9/225-t10/225); + const double t35 = (-t1/900-2.0/225.0*t3-7.0/900.0*t4+7.0/900.0*t5 + +t6/900+7.0/900.0*t12 + -t7/900+2.0/225.0*t8-7.0/900.0*t13+t10/900); + const double t36 = (t2/900-7.0/900.0*t3-2.0/225.0*t4+2.0/225.0*t5 + +7.0/900.0*t12+t7/900+ + 7.0/900.0*t8-7.0/900.0*t13-t9/900-t10/900); + const double t37 = (-t3/225-t4/225+t5/225+t12/225+t8/225-t13/225); + const double t38 = (-14.0/225.0*t1+8.0/225.0*t2-8.0/225.0*t3 + -2.0/225.0*t4+2.0/225.0*t5+ + 14.0/225.0*t6-2.0/75.0*t12-2.0/75.0*t7 + +8.0/225.0*t8+2.0/75.0*t13-8.0/225.0*t9+ + 2.0/75.0*t10); + const double t39 = (2.0/225.0*t1-2.0/225.0*t2+2.0/225.0*t3 + +2.0/225.0*t4-2.0/225.0*t5 + -2.0/225.0*t6-2.0/225.0*t8+2.0/225.0*t9); + const double t40 = (-8.0/225.0*t1+4.0/225.0*t2-4.0/225.0*t3 + +8.0/225.0*t6-4.0/225.0*t12 + -4.0/225.0*t7+4.0/225.0*t8+4.0/225.0*t13 + -4.0/225.0*t9+4.0/225.0*t10); + const double t41 = (-8.0/225.0*t1+14.0/225.0*t2-2.0/225.0*t3 + -8.0/225.0*t4+8.0/225.0*t5+ + 8.0/225.0*t6-2.0/75.0*t12+2.0/75.0*t7 + +2.0/225.0*t8+2.0/75.0*t13-14.0/225.0*t9 + -2.0/75.0*t10); + const double t42 = (-4.0/225.0*t1+8.0/225.0*t2-4.0/225.0*t4 + +4.0/225.0*t5+4.0/225.0*t6 + -4.0/225.0*t12+4.0/225.0*t7+4.0/225.0*t13 + -8.0/225.0*t9-4.0/225.0*t10); + const double t43 = (-2.0/225.0*t1+8.0/225.0*t2-8.0/225.0*t3 + -14.0/225.0*t4+14.0/225.0*t5 + +2.0/225.0*t6+2.0/75.0*t12+2.0/75.0*t7 + +8.0/225.0*t8-2.0/75.0*t13-8.0/225.0*t9 + -2.0/75.0*t10); + const double t44 = (4.0/225.0*t2-4.0/225.0*t3-8.0/225.0*t4 + +8.0/225.0*t5+4.0/225.0*t12+ + 4.0/225.0*t7+4.0/225.0*t8-4.0/225.0*t13 + -4.0/225.0*t9-4.0/225.0*t10); + const double t45 = (-8.0/225.0*t1+2.0/225.0*t2-14.0/225.0*t3 + -8.0/225.0*t4+8.0/225.0*t5+ + 8.0/225.0*t6+2.0/75.0*t12-2.0/75.0*t7 + +14.0/225.0*t8-2.0/75.0*t13-2.0/225.0*t9+ + 2.0/75.0*t10); + const double t46 = (-4.0/225.0*t1-8.0/225.0*t3-4.0/225.0*t4 + +4.0/225.0*t5+4.0/225.0*t6+ + 4.0/225.0*t12-4.0/225.0*t7+8.0/225.0*t8 + -4.0/225.0*t13+4.0/225.0*t10); + + local_mass_matrix(0,0) = (-7.0/450.0*t1+t2/450-7.0/450.0*t3 + -t4/450+t5/450+7.0/450.0*t6-t7/75 + +7.0/450.0*t8-t9/450+t10/75); + local_mass_matrix(0,1) = (t14); + local_mass_matrix(0,2) = (t15); + local_mass_matrix(0,3) = (t16); + local_mass_matrix(0,4) = (t17); + local_mass_matrix(0,5) = (t18); + local_mass_matrix(0,6) = (t19); + local_mass_matrix(0,7) = (t20); + local_mass_matrix(0,8) = (t21); + local_mass_matrix(1,0) = (t14); + local_mass_matrix(1,1) = (-7.0/450.0*t1+7.0/450.0*t2-t3/450 + -t4/450+t5/450+7.0/450.0*t6- + t12/75+t8/450+t13/75-7.0/450.0*t9); + local_mass_matrix(1,2) = (t23); + local_mass_matrix(1,3) = (t15); + local_mass_matrix(1,4) = (t24); + local_mass_matrix(1,5) = (t25); + local_mass_matrix(1,6) = (t26); + local_mass_matrix(1,7) = (t18); + local_mass_matrix(1,8) = (t27); + local_mass_matrix(2,0) = (t15); + local_mass_matrix(2,1) = (t23); + local_mass_matrix(2,2) = (-t1/450+7.0/450.0*t2-t3/450-7.0/450.0*t4 + +7.0/450.0*t5+t6/450+t7/75 + +t8/450-7.0/450.0*t9-t10/75); + local_mass_matrix(2,3) = (t29); + local_mass_matrix(2,4) = (t26); + local_mass_matrix(2,5) = (t30); + local_mass_matrix(2,6) = (t31); + local_mass_matrix(2,7) = (t32); + local_mass_matrix(2,8) = (t33); + local_mass_matrix(3,0) = (t16); + local_mass_matrix(3,1) = (t15); + local_mass_matrix(3,2) = (t29); + local_mass_matrix(3,3) = (-t1/450+t2/450-7.0/450.0*t3-7.0/450.0*t4 + +7.0/450.0*t5+t6/450+ + t12/75+7.0/450.0*t8-t13/75-t9/450); + local_mass_matrix(3,4) = (t19); + local_mass_matrix(3,5) = (t32); + local_mass_matrix(3,6) = (t35); + local_mass_matrix(3,7) = (t36); + local_mass_matrix(3,8) = (t37); + local_mass_matrix(4,0) = (t17); + local_mass_matrix(4,1) = (t24); + local_mass_matrix(4,2) = (t26); + local_mass_matrix(4,3) = (t19); + local_mass_matrix(4,4) = (t38); + local_mass_matrix(4,5) = (t27); + local_mass_matrix(4,6) = (t39); + local_mass_matrix(4,7) = (t21); + local_mass_matrix(4,8) = (t40); + local_mass_matrix(5,0) = (t18); + local_mass_matrix(5,1) = (t25); + local_mass_matrix(5,2) = (t30); + local_mass_matrix(5,3) = (t32); + local_mass_matrix(5,4) = (t27); + local_mass_matrix(5,5) = (t41); + local_mass_matrix(5,6) = (t33); + local_mass_matrix(5,7) = (t39); + local_mass_matrix(5,8) = (t42); + local_mass_matrix(6,0) = (t19); + local_mass_matrix(6,1) = (t26); + local_mass_matrix(6,2) = (t31); + local_mass_matrix(6,3) = (t35); + local_mass_matrix(6,4) = (t39); + local_mass_matrix(6,5) = (t33); + local_mass_matrix(6,6) = (t43); + local_mass_matrix(6,7) = (t37); + local_mass_matrix(6,8) = (t44); + local_mass_matrix(7,0) = (t20); + local_mass_matrix(7,1) = (t18); + local_mass_matrix(7,2) = (t32); + local_mass_matrix(7,3) = (t36); + local_mass_matrix(7,4) = (t21); + local_mass_matrix(7,5) = (t39); + local_mass_matrix(7,6) = (t37); + local_mass_matrix(7,7) = (t45); + local_mass_matrix(7,8) = (t46); + local_mass_matrix(8,0) = (t21); + local_mass_matrix(8,1) = (t27); + local_mass_matrix(8,2) = (t33); + local_mass_matrix(8,3) = (t37); + local_mass_matrix(8,4) = (t40); + local_mass_matrix(8,5) = (t42); + local_mass_matrix(8,6) = (t44); + local_mass_matrix(8,7) = (t46); + local_mass_matrix(8,8) = (-32.0/225.0*t1+32.0/225.0*t2-32.0/225.0*t3 + -32.0/225.0*t4+32.0/225.0*t5+32.0/225.0*t6 + +32.0/225.0*t8-32.0/225.0*t9); +}; + + + template <> void FEQuadraticSub<2>::get_ansatz_points (const typename DoFHandler<2>::cell_iterator &cell, const Boundary<2>&, @@ -1036,19 +1335,6 @@ void FEQuadraticSub::fill_fe_values (const DoFHandler::cell_iterator & - -template -void FEQuadraticSub::get_local_mass_matrix (const DoFHandler::cell_iterator &, - const Boundary &, - dFMatrix &) const { - Assert (false, ExcNotImplemented()); -}; - - - - - - #if deal_II_dimension == 1 template <> -- 2.39.5