From f124a131c3381d7c5b142b6de41a5109553daa8d Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 26 Jun 1998 12:24:15 +0000 Subject: [PATCH] Split fe_lib.cc into linear and quadratic elements. git-svn-id: https://svn.dealii.org/trunk@413 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_lib.linear.cc | 924 +---------------- deal.II/deal.II/source/fe/fe_lib.quadratic.cc | 929 ++++++++++++++++++ 2 files changed, 934 insertions(+), 919 deletions(-) create mode 100644 deal.II/deal.II/source/fe/fe_lib.quadratic.cc 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 11ef7c4980..a58dd9cd12 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -7,6 +7,8 @@ #include + + #if deal_II_dimension == 1 template <> @@ -270,6 +272,7 @@ FELinear<2>::FELinear () : template <> +inline double FELinear<2>::shape_value (const unsigned int i, const Point<2>& p) const @@ -559,7 +562,8 @@ void FELinear::fill_fe_values (const DoFHandler::cell_iterator &cell, // x_l(xi_l) = sum_j p_j N_j(xi_l) for (unsigned int j=0; j::vertices_per_cell; ++j) for (unsigned int l=0; l::shape_value(j, unit_points[l]); }; @@ -641,924 +645,6 @@ void FELinear::get_face_ansatz_points (const typename DoFHandler::face - - - -#if deal_II_dimension == 1 - -template <> -FEQuadraticSub<1>::FEQuadraticSub () : - FiniteElement<1> (1, 1) {}; - - - -template <> -void FEQuadraticSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - const Boundary<1> &boundary) const { - // simply pass down - FiniteElement<1>::fill_fe_values (cell, unit_points, - jacobians, compute_jacobians, - ansatz_points, compute_ansatz_points, - q_points, compute_q_points, boundary); -}; - - - -template <> -double -FEQuadraticSub<1>::shape_value(const unsigned int i, - const Point<1> &p) const -{ - Assert((i -inline -Point<1> -FEQuadraticSub<1>::shape_grad(const unsigned int i, - const Point<1> &p) const -{ - Assert((i(-3+4*xi); - case 1: return Point<1>(4*xi-1); - case 2: return Point<1>(4-8*xi); - } - return Point<1>(); -}; - - - -template <> -void FEQuadraticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell, - const Boundary<1> &boundary, - vector > &ansatz_points) const { - FiniteElement<1>::get_ansatz_points (cell, boundary, ansatz_points); -}; - - - -template <> -void FEQuadraticSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &, - const Boundary<1> &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, - const Boundary<1> &, - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, - const unsigned int , - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const Boundary<1> &, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const unsigned int, - const vector > &, - vector > &) const { - 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 - - -#if deal_II_dimension == 2 - -template <> -FEQuadraticSub<2>::FEQuadraticSub () : - FiniteElement<2> (1, 1, 1) -{ - interface_constraints(0,2) = 1.0; - interface_constraints(1,0) = 3./8.; - interface_constraints(1,1) = -1./8.; - interface_constraints(1,2) = 3./4.; - interface_constraints(2,0) = -1./8.; - interface_constraints(2,1) = 3./8.; - interface_constraints(2,2) = 3./4.; - - // still implement restriction - // and prolongation - Assert (false, ExcNotImplemented()); -}; - - -template <> -double -FEQuadraticSub<2>::shape_value (const unsigned int i, - const Point<2> &p) const -{ - Assert (i -Point<2> -FEQuadraticSub<2>::shape_grad (const unsigned int i, - const Point<2> &p) const -{ - Assert (i(-(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>&, - vector > &ansatz_points) const { - Assert (ansatz_points.size() == total_dofs, - ExcWrongFieldDimension (ansatz_points.size(), total_dofs)); - - for (unsigned int vertex=0; vertex<4; ++vertex) - ansatz_points[vertex] = cell->vertex(vertex); - - // for the bilinear mapping, the centers - // of the face on the unit cell are mapped - // to the mean coordinates of the vertices - for (unsigned int line=0; line<4; ++line) - ansatz_points[4+line] = (cell->line(line)->vertex(0) + - cell->line(line)->vertex(1)) / 2; - // same for the center of the square: - // since all four linear basis functions - // take on the value 1/4 at the center, - // the center is mapped to the mean - // coordinates of the four vertices - ansatz_points[8] = (ansatz_points[0] + - ansatz_points[1] + - ansatz_points[2] + - ansatz_points[3]) / 4; -}; - - - -template <> -void FEQuadraticSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::face_iterator &face, - const Boundary<2> &, - vector > &ansatz_points) const { - Assert (ansatz_points.size() == dofs_per_face, - ExcWrongFieldDimension (ansatz_points.size(), dofs_per_face)); - - for (unsigned int vertex=0; vertex<2; ++vertex) - ansatz_points[vertex] = face->vertex(vertex); - ansatz_points[2] = (ansatz_points[0] + ansatz_points[1]) / 1; -}; - - - -template <> -void FEQuadraticSub<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face, - const Boundary<2> &, - const vector > &unit_points, - vector &face_jacobians) const { - // more or less copied from the linear - // finite element - Assert (unit_points.size() == face_jacobians.size(), - ExcWrongFieldDimension (unit_points.size(), face_jacobians.size())); - - // a linear mapping for a single line - // produces particularly simple - // expressions for the jacobi - // determinant :-) - const double h = sqrt((face->vertex(1) - face->vertex(0)).square()); - fill_n (face_jacobians.begin(), - unit_points.size(), - h); -}; - - - -template <> -void FEQuadraticSub<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face, - const unsigned int , - const vector > &unit_points, - vector &face_jacobians) const { - // more or less copied from the linear - // finite element - Assert (unit_points.size() == face_jacobians.size(), - ExcWrongFieldDimension (unit_points.size(), face_jacobians.size())); - Assert (face->at_boundary() == false, - ExcBoundaryFaceUsed ()); - - // a linear mapping for a single line - // produces particularly simple - // expressions for the jacobi - // determinant :-) - const double h = sqrt((face->vertex(1) - face->vertex(0)).square()); - fill_n (face_jacobians.begin(), - unit_points.size(), - h/2); -}; - - - -template <> -void FEQuadraticSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell, - const unsigned int face_no, - const Boundary<2> &, - const vector > &unit_points, - vector > &normal_vectors) const { - // more or less copied from the linear - // finite element - Assert (unit_points.size() == normal_vectors.size(), - ExcWrongFieldDimension (unit_points.size(), normal_vectors.size())); - - const DoFHandler<2>::face_iterator face = cell->face(face_no); - // compute direction of line - const Point<2> line_direction = (face->vertex(1) - face->vertex(0)); - // rotate to the right by 90 degrees - const Point<2> normal_direction(line_direction(1), - -line_direction(0)); - - if (face_no <= 1) - // for sides 0 and 1: return the correctly - // scaled vector - fill (normal_vectors.begin(), normal_vectors.end(), - normal_direction / sqrt(normal_direction.square())); - else - // for sides 2 and 3: scale and invert - // vector - fill (normal_vectors.begin(), normal_vectors.end(), - normal_direction / (-sqrt(normal_direction.square()))); -}; - - - -template <> -void FEQuadraticSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell, - const unsigned int face_no, - const unsigned int, - const vector > &unit_points, - vector > &normal_vectors) const { - // more or less copied from the linear - // finite element - // note, that in 2D the normal vectors to the - // subface have the same direction as that - // for the face - Assert (unit_points.size() == normal_vectors.size(), - ExcWrongFieldDimension (unit_points.size(), normal_vectors.size())); - Assert (cell->face(face_no)->at_boundary() == false, - ExcBoundaryFaceUsed ()); - - const DoFHandler<2>::face_iterator face = cell->face(face_no); - // compute direction of line - const Point<2> line_direction = (face->vertex(1) - face->vertex(0)); - // rotate to the right by 90 degrees - const Point<2> normal_direction(line_direction(1), - -line_direction(0)); - - if (face_no <= 1) - // for sides 0 and 1: return the correctly - // scaled vector - fill (normal_vectors.begin(), normal_vectors.end(), - normal_direction / sqrt(normal_direction.square())); - else - // for sides 2 and 3: scale and invert - // vector - fill (normal_vectors.begin(), normal_vectors.end(), - normal_direction / (-sqrt(normal_direction.square()))); -}; - -#endif - - - - - -template -void FEQuadraticSub::fill_fe_values (const DoFHandler::cell_iterator &, - const vector > &unit_points, - vector &jacobians, - const bool, - vector > &ansatz_points, - const bool, - vector > &q_points, - const bool, - const Boundary &) const { - Assert (jacobians.size() == unit_points.size(), - ExcWrongFieldDimension(jacobians.size(), unit_points.size())); - Assert (q_points.size() == unit_points.size(), - ExcWrongFieldDimension(q_points.size(), unit_points.size())); - Assert (ansatz_points.size() == total_dofs, - ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); - - Assert (false, ExcNotImplemented()); -}; - - - - - - -#if deal_II_dimension == 1 - -template <> -FECubic<1>::FECubic () : - FiniteElement<1> (1, 2) {}; - - - -template <> -void FECubic<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, - const vector > &unit_points, - vector &jacobians, - const bool compute_jacobians, - vector > &ansatz_points, - const bool compute_ansatz_points, - vector > &q_points, - const bool compute_q_points, - const Boundary<1> &boundary) const { - // simply pass down - FiniteElement<1>::fill_fe_values (cell, unit_points, - jacobians, compute_jacobians, - ansatz_points, compute_ansatz_points, - q_points, compute_q_points, boundary); -}; - - - -template <> -void FECubic<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell, - const Boundary<1> &boundary, - vector > &ansatz_points) const { - FiniteElement<1>::get_ansatz_points (cell, boundary, ansatz_points); -}; - - - -template <> -void FECubic<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &, - const Boundary<1> &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FECubic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, - const Boundary<1> &, - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FECubic<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, - const unsigned int , - const vector > &, - vector &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FECubic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const Boundary<1> &, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - - - -template <> -void FECubic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, - const unsigned int, - const unsigned int, - const vector > &, - vector > &) const { - Assert (false, ExcInternalError()); -}; - -#endif - - - -#if deal_II_dimension == 2 - -template <> -FECubic<2>::FECubic () : - FiniteElement<2> (1, 2, 4) {}; - -#endif - - - -template -double -FECubic::shape_value (const unsigned int i, - const Point &) const -{ - Assert (i::ExcInvalidIndex(i)); - Assert (false, ExcNotImplemented()); - return 0.; -}; - - - -template -Point -FECubic::shape_grad (const unsigned int i, - const Point &) const -{ - Assert (i::ExcInvalidIndex(i)); - Assert (false, ExcNotImplemented()); - return Point (); -}; - - - -template -void FECubic::fill_fe_values (const DoFHandler::cell_iterator &, - const vector > &unit_points, - vector &jacobians, - const bool, - vector > &ansatz_points, - const bool, - vector > &q_points, - const bool, - const Boundary &) const { - Assert (jacobians.size() == unit_points.size(), - ExcWrongFieldDimension(jacobians.size(), unit_points.size())); - Assert (q_points.size() == unit_points.size(), - ExcWrongFieldDimension(q_points.size(), unit_points.size())); - Assert (ansatz_points.size() == total_dofs, - ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); - - Assert (false, ExcNotImplemented()); -}; - - - -template -void FECubic::get_ansatz_points (const typename DoFHandler::cell_iterator &, - const Boundary &, - vector > &) const { - Assert (false, ExcNotImplemented()); -}; - - - -template -void FECubic::get_face_ansatz_points (const typename DoFHandler::face_iterator &, - const Boundary &, - vector > &) const { - Assert (false, ExcNotImplemented()); -}; - - - -template -void FECubic::get_face_jacobians (const DoFHandler::face_iterator &, - const Boundary &, - const vector > &, - vector &) const { - Assert (false, ExcNotImplemented()); -}; - - - -template -void FECubic::get_subface_jacobians (const DoFHandler::face_iterator &face, - const unsigned int , - const vector > &, - vector &) const { - Assert (face->at_boundary() == false, - ExcBoundaryFaceUsed ()); - - Assert (false, ExcNotImplemented()); -}; - - - -template -void FECubic::get_normal_vectors (const DoFHandler::cell_iterator &, - const unsigned int, - const Boundary &, - const vector > &, - vector > &) const { - Assert (false, ExcNotImplemented()); -}; - - - -template -void FECubic::get_normal_vectors (const DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const unsigned int , - const vector > &, - vector > &) const { - Assert (cell->face(face_no)->at_boundary() == false, - ExcBoundaryFaceUsed ()); - - Assert (false, ExcNotImplemented()); -}; - - - -template -void FECubic::get_local_mass_matrix (const DoFHandler::cell_iterator &, - const Boundary &, - dFMatrix &) const { - Assert (false, ExcNotImplemented()); -}; - - - - // explicit instantiations template class FELinear; -template class FEQuadraticSub; -template class FECubic; - diff --git a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc new file mode 100644 index 0000000000..49ec2a0745 --- /dev/null +++ b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc @@ -0,0 +1,929 @@ +/* $Id$ */ + +#include +#include +#include +#include +#include + + + + + + +#if deal_II_dimension == 1 + +template <> +FEQuadraticSub<1>::FEQuadraticSub () : + FiniteElement<1> (1, 1) {}; + + + +template <> +void FEQuadraticSub<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, + const vector > &unit_points, + vector &jacobians, + const bool compute_jacobians, + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points, + const Boundary<1> &boundary) const { + // simply pass down + FiniteElement<1>::fill_fe_values (cell, unit_points, + jacobians, compute_jacobians, + ansatz_points, compute_ansatz_points, + q_points, compute_q_points, boundary); +}; + + + +template <> +double +FEQuadraticSub<1>::shape_value(const unsigned int i, + const Point<1> &p) const +{ + Assert((i +inline +Point<1> +FEQuadraticSub<1>::shape_grad(const unsigned int i, + const Point<1> &p) const +{ + Assert((i(-3+4*xi); + case 1: return Point<1>(4*xi-1); + case 2: return Point<1>(4-8*xi); + } + return Point<1>(); +}; + + + +template <> +void FEQuadraticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell, + const Boundary<1> &boundary, + vector > &ansatz_points) const { + FiniteElement<1>::get_ansatz_points (cell, boundary, ansatz_points); +}; + + + +template <> +void FEQuadraticSub<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &, + const Boundary<1> &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FEQuadraticSub<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, + const Boundary<1> &, + const vector > &, + vector &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FEQuadraticSub<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, + const unsigned int , + const vector > &, + vector &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const Boundary<1> &, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FEQuadraticSub<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const unsigned int, + const vector > &, + vector > &) const { + 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 + + +#if deal_II_dimension == 2 + +template <> +FEQuadraticSub<2>::FEQuadraticSub () : + FiniteElement<2> (1, 1, 1) +{ + interface_constraints(0,2) = 1.0; + interface_constraints(1,0) = 3./8.; + interface_constraints(1,1) = -1./8.; + interface_constraints(1,2) = 3./4.; + interface_constraints(2,0) = -1./8.; + interface_constraints(2,1) = 3./8.; + interface_constraints(2,2) = 3./4.; + + // still implement restriction + // and prolongation + Assert (false, ExcNotImplemented()); +}; + + +template <> +double +FEQuadraticSub<2>::shape_value (const unsigned int i, + const Point<2> &p) const +{ + Assert (i +Point<2> +FEQuadraticSub<2>::shape_grad (const unsigned int i, + const Point<2> &p) const +{ + Assert (i(-(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>&, + vector > &ansatz_points) const { + Assert (ansatz_points.size() == total_dofs, + ExcWrongFieldDimension (ansatz_points.size(), total_dofs)); + + for (unsigned int vertex=0; vertex<4; ++vertex) + ansatz_points[vertex] = cell->vertex(vertex); + + // for the bilinear mapping, the centers + // of the face on the unit cell are mapped + // to the mean coordinates of the vertices + for (unsigned int line=0; line<4; ++line) + ansatz_points[4+line] = (cell->line(line)->vertex(0) + + cell->line(line)->vertex(1)) / 2; + // same for the center of the square: + // since all four linear basis functions + // take on the value 1/4 at the center, + // the center is mapped to the mean + // coordinates of the four vertices + ansatz_points[8] = (ansatz_points[0] + + ansatz_points[1] + + ansatz_points[2] + + ansatz_points[3]) / 4; +}; + + + +template <> +void FEQuadraticSub<2>::get_face_ansatz_points (const typename DoFHandler<2>::face_iterator &face, + const Boundary<2> &, + vector > &ansatz_points) const { + Assert (ansatz_points.size() == dofs_per_face, + ExcWrongFieldDimension (ansatz_points.size(), dofs_per_face)); + + for (unsigned int vertex=0; vertex<2; ++vertex) + ansatz_points[vertex] = face->vertex(vertex); + ansatz_points[2] = (ansatz_points[0] + ansatz_points[1]) / 1; +}; + + + +template <> +void FEQuadraticSub<2>::get_face_jacobians (const DoFHandler<2>::face_iterator &face, + const Boundary<2> &, + const vector > &unit_points, + vector &face_jacobians) const { + // more or less copied from the linear + // finite element + Assert (unit_points.size() == face_jacobians.size(), + ExcWrongFieldDimension (unit_points.size(), face_jacobians.size())); + + // a linear mapping for a single line + // produces particularly simple + // expressions for the jacobi + // determinant :-) + const double h = sqrt((face->vertex(1) - face->vertex(0)).square()); + fill_n (face_jacobians.begin(), + unit_points.size(), + h); +}; + + + +template <> +void FEQuadraticSub<2>::get_subface_jacobians (const DoFHandler<2>::face_iterator &face, + const unsigned int , + const vector > &unit_points, + vector &face_jacobians) const { + // more or less copied from the linear + // finite element + Assert (unit_points.size() == face_jacobians.size(), + ExcWrongFieldDimension (unit_points.size(), face_jacobians.size())); + Assert (face->at_boundary() == false, + ExcBoundaryFaceUsed ()); + + // a linear mapping for a single line + // produces particularly simple + // expressions for the jacobi + // determinant :-) + const double h = sqrt((face->vertex(1) - face->vertex(0)).square()); + fill_n (face_jacobians.begin(), + unit_points.size(), + h/2); +}; + + + +template <> +void FEQuadraticSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell, + const unsigned int face_no, + const Boundary<2> &, + const vector > &unit_points, + vector > &normal_vectors) const { + // more or less copied from the linear + // finite element + Assert (unit_points.size() == normal_vectors.size(), + ExcWrongFieldDimension (unit_points.size(), normal_vectors.size())); + + const DoFHandler<2>::face_iterator face = cell->face(face_no); + // compute direction of line + const Point<2> line_direction = (face->vertex(1) - face->vertex(0)); + // rotate to the right by 90 degrees + const Point<2> normal_direction(line_direction(1), + -line_direction(0)); + + if (face_no <= 1) + // for sides 0 and 1: return the correctly + // scaled vector + fill (normal_vectors.begin(), normal_vectors.end(), + normal_direction / sqrt(normal_direction.square())); + else + // for sides 2 and 3: scale and invert + // vector + fill (normal_vectors.begin(), normal_vectors.end(), + normal_direction / (-sqrt(normal_direction.square()))); +}; + + + +template <> +void FEQuadraticSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator &cell, + const unsigned int face_no, + const unsigned int, + const vector > &unit_points, + vector > &normal_vectors) const { + // more or less copied from the linear + // finite element + // note, that in 2D the normal vectors to the + // subface have the same direction as that + // for the face + Assert (unit_points.size() == normal_vectors.size(), + ExcWrongFieldDimension (unit_points.size(), normal_vectors.size())); + Assert (cell->face(face_no)->at_boundary() == false, + ExcBoundaryFaceUsed ()); + + const DoFHandler<2>::face_iterator face = cell->face(face_no); + // compute direction of line + const Point<2> line_direction = (face->vertex(1) - face->vertex(0)); + // rotate to the right by 90 degrees + const Point<2> normal_direction(line_direction(1), + -line_direction(0)); + + if (face_no <= 1) + // for sides 0 and 1: return the correctly + // scaled vector + fill (normal_vectors.begin(), normal_vectors.end(), + normal_direction / sqrt(normal_direction.square())); + else + // for sides 2 and 3: scale and invert + // vector + fill (normal_vectors.begin(), normal_vectors.end(), + normal_direction / (-sqrt(normal_direction.square()))); +}; + +#endif + + + + + +template +void FEQuadraticSub::fill_fe_values (const DoFHandler::cell_iterator &, + const vector > &unit_points, + vector &jacobians, + const bool, + vector > &ansatz_points, + const bool, + vector > &q_points, + const bool, + const Boundary &) const { + Assert (jacobians.size() == unit_points.size(), + ExcWrongFieldDimension(jacobians.size(), unit_points.size())); + Assert (q_points.size() == unit_points.size(), + ExcWrongFieldDimension(q_points.size(), unit_points.size())); + Assert (ansatz_points.size() == total_dofs, + ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); + + Assert (false, ExcNotImplemented()); +}; + + + + + + +#if deal_II_dimension == 1 + +template <> +FECubic<1>::FECubic () : + FiniteElement<1> (1, 2) {}; + + + +template <> +void FECubic<1>::fill_fe_values (const DoFHandler<1>::cell_iterator &cell, + const vector > &unit_points, + vector &jacobians, + const bool compute_jacobians, + vector > &ansatz_points, + const bool compute_ansatz_points, + vector > &q_points, + const bool compute_q_points, + const Boundary<1> &boundary) const { + // simply pass down + FiniteElement<1>::fill_fe_values (cell, unit_points, + jacobians, compute_jacobians, + ansatz_points, compute_ansatz_points, + q_points, compute_q_points, boundary); +}; + + + +template <> +void FECubic<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell, + const Boundary<1> &boundary, + vector > &ansatz_points) const { + FiniteElement<1>::get_ansatz_points (cell, boundary, ansatz_points); +}; + + + +template <> +void FECubic<1>::get_face_ansatz_points (const typename DoFHandler<1>::face_iterator &, + const Boundary<1> &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FECubic<1>::get_face_jacobians (const DoFHandler<1>::face_iterator &, + const Boundary<1> &, + const vector > &, + vector &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FECubic<1>::get_subface_jacobians (const DoFHandler<1>::face_iterator &, + const unsigned int , + const vector > &, + vector &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FECubic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const Boundary<1> &, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + + + +template <> +void FECubic<1>::get_normal_vectors (const DoFHandler<1>::cell_iterator &, + const unsigned int, + const unsigned int, + const vector > &, + vector > &) const { + Assert (false, ExcInternalError()); +}; + +#endif + + + +#if deal_II_dimension == 2 + +template <> +FECubic<2>::FECubic () : + FiniteElement<2> (1, 2, 4) {}; + +#endif + + + +template +double +FECubic::shape_value (const unsigned int i, + const Point &) const +{ + Assert (i::ExcInvalidIndex(i)); + Assert (false, ExcNotImplemented()); + return 0.; +}; + + + +template +Point +FECubic::shape_grad (const unsigned int i, + const Point &) const +{ + Assert (i::ExcInvalidIndex(i)); + Assert (false, ExcNotImplemented()); + return Point (); +}; + + + +template +void FECubic::fill_fe_values (const DoFHandler::cell_iterator &, + const vector > &unit_points, + vector &jacobians, + const bool, + vector > &ansatz_points, + const bool, + vector > &q_points, + const bool, + const Boundary &) const { + Assert (jacobians.size() == unit_points.size(), + ExcWrongFieldDimension(jacobians.size(), unit_points.size())); + Assert (q_points.size() == unit_points.size(), + ExcWrongFieldDimension(q_points.size(), unit_points.size())); + Assert (ansatz_points.size() == total_dofs, + ExcWrongFieldDimension(ansatz_points.size(), total_dofs)); + + Assert (false, ExcNotImplemented()); +}; + + + +template +void FECubic::get_ansatz_points (const typename DoFHandler::cell_iterator &, + const Boundary &, + vector > &) const { + Assert (false, ExcNotImplemented()); +}; + + + +template +void FECubic::get_face_ansatz_points (const typename DoFHandler::face_iterator &, + const Boundary &, + vector > &) const { + Assert (false, ExcNotImplemented()); +}; + + + +template +void FECubic::get_face_jacobians (const DoFHandler::face_iterator &, + const Boundary &, + const vector > &, + vector &) const { + Assert (false, ExcNotImplemented()); +}; + + + +template +void FECubic::get_subface_jacobians (const DoFHandler::face_iterator &face, + const unsigned int , + const vector > &, + vector &) const { + Assert (face->at_boundary() == false, + ExcBoundaryFaceUsed ()); + + Assert (false, ExcNotImplemented()); +}; + + + +template +void FECubic::get_normal_vectors (const DoFHandler::cell_iterator &, + const unsigned int, + const Boundary &, + const vector > &, + vector > &) const { + Assert (false, ExcNotImplemented()); +}; + + + +template +void FECubic::get_normal_vectors (const DoFHandler::cell_iterator &cell, + const unsigned int face_no, + const unsigned int , + const vector > &, + vector > &) const { + Assert (cell->face(face_no)->at_boundary() == false, + ExcBoundaryFaceUsed ()); + + Assert (false, ExcNotImplemented()); +}; + + + +template +void FECubic::get_local_mass_matrix (const DoFHandler::cell_iterator &, + const Boundary &, + dFMatrix &) const { + Assert (false, ExcNotImplemented()); +}; + + + + +// explicit instantiations + +template class FEQuadraticSub; + -- 2.39.5