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
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;
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 >t<phi), while the phi[i]
+ are the basis functions of the biquadratic element.
+
+ x_real := sum(x[i]*tphi[i], i=0..3);
+ y_real := sum(y[i]*tphi[i], i=0..3);
+ tphi[0] := (1-xi)*(1-eta);
+ tphi[1] := xi*(1-eta);
+ tphi[2] := xi*eta;
+ tphi[3] := (1-xi)*eta;
+ detJ := diff(x_real,xi)*diff(y_real,eta) - diff(x_real,eta)*diff(y_real,xi);
+
+ phi[0] := (1-xi)*( 2*xi-1) * (1-eta)*( 2*eta-1);
+ phi[1] := xi *(-2*xi+1) * (1-eta)*( 2*eta-1);
+ phi[2] := xi *(-2*xi+1) * eta *(-2*eta+1);
+ phi[3] := (1-xi)*( 2*xi-1) * eta *(-2*eta+1);
+ phi[4] := 4 * (1-xi)*xi * (1-eta)*(1-2*eta);
+ phi[5] := 4 * xi *(-1+2*xi) * (1-eta)*eta;
+ phi[6] := 4 * (1-xi)*xi * eta *(-1+2*eta);
+ phi[7] := 4 * (1-xi)*(1-2*xi) * (1-eta)*eta;
+ phi[8] := 16 * xi*(1-xi) * eta*(1-eta);
+ m := proc (i,j) int( int(phi[i]*phi[j]*detJ, xi=0..1), eta=0..1); end;
+
+ M := array(0..8,0..8);
+ for i from 0 to 8 do
+ for j from 0 to 8 do
+ M[i,j] := m(i,j);
+ od;
+ od;
+
+ readlib(C);
+ C(M, optimized);
+*/
+
+ const double x[4] = { cell->vertex(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>&,
-
-template <int dim>
-void FEQuadraticSub<dim>::get_local_mass_matrix (const DoFHandler<dim>::cell_iterator &,
- const Boundary<dim> &,
- dFMatrix &) const {
- Assert (false, ExcNotImplemented());
-};
-
-
-
-
-
-
#if deal_II_dimension == 1
template <>