From 8ae4e4915f295e33602d90bf10c361d6bb51e2e2 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 26 Jun 1998 16:20:19 +0000 Subject: [PATCH] Further on the road with implementation of quadratic elements. git-svn-id: https://svn.dealii.org/trunk@414 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/Todo | 3 + deal.II/deal.II/include/fe/fe.h | 14 +- deal.II/deal.II/include/fe/fe_lib.lagrange.h | 116 +-- deal.II/deal.II/source/fe/fe_lib.quadratic.cc | 754 ++++++++++++------ 4 files changed, 551 insertions(+), 336 deletions(-) diff --git a/deal.II/deal.II/Todo b/deal.II/deal.II/Todo index 8b1386520b..0645b9a276 100644 --- a/deal.II/deal.II/Todo +++ b/deal.II/deal.II/Todo @@ -118,6 +118,9 @@ Fully implement the POVRAY format, i.e. use textures, a better angle of view, etc. The present implementation is only a rudimentary hack. +Review the restruction matrices. I'm not really sure about their + meaning and how they are defined, so they may be wrong for linear + elements and they are not implemented at all for quadratic ones. diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 11df68192d..fa843db075 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -306,11 +306,23 @@ struct FiniteElementBase : public FiniteElementData { * are for the unrefined cell's degrees of * freedom. Thus, if #u0# is the vector * of values of degrees of freedom on the - * coarse cell, * #prolongation[i]*u0# + * coarse cell, #prolongation[i]*u0# * yields the vector of values of the * degrees of freedom on the #i#th child * cell. * + * On the other hand, for finite elements + * with embedded spaces, the basis function + * phi0[i] on the coarse grid can be + * expressed by + * \sum_c \sum_j p^c_{ji) phi1[j] + * where the sum over c runs over the child + * cells and phi1[j] is the j_th basis + * function on the c_th child cell. Note + * that we need here the transpose of the + * matrix p^c (p^c is returned by this + * function with parameter c). + * * Upon assembling the transfer matrix * between cells using this matrix array, * zero elements in the prolongation diff --git a/deal.II/deal.II/include/fe/fe_lib.lagrange.h b/deal.II/deal.II/include/fe/fe_lib.lagrange.h index 1b3f403463..8d54ce4c61 100644 --- a/deal.II/deal.II/include/fe/fe_lib.lagrange.h +++ b/deal.II/deal.II/include/fe/fe_lib.lagrange.h @@ -269,118 +269,32 @@ class FEQuadraticSub : public FiniteElement { virtual void get_local_mass_matrix (const DoFHandler::cell_iterator &cell, const Boundary &boundary, dFMatrix &local_mass_matrix) const; -}; - - - - -/** - * Define a (bi-, tri-, etc)cubic finite element in #dim# space dimensions. - * In one space dimension, a linear (subparametric) mapping from the unit cell - * to the real cell is implemented. - */ -template -class FECubic : public FiniteElement { - public: - /** - * Constructor - */ - FECubic (); - /** + private: + /** * Return the value of the #i#th shape * function at point #p# on the unit cell. + * Here, the (bi-)linear basis functions + * are meant, which are used for the + * computation of the transformation from + * unit cell to real space cell. */ - virtual double shape_value(const unsigned int i, - const Point& p) const; + double linear_shape_value(const unsigned int i, + const Point& p) const; /** * Return the gradient of the #i#th shape * function at point #p# on the unit cell. + * Here, the (bi-)linear basis functions + * are meant, which are used for the + * computation of the transformation from + * unit cell to real space cell. */ - virtual Point shape_grad(const unsigned int i, - const Point& p) const; - - /** - * Refer to the base class for detailed - * information on this function. - * - * For one dimensional elements, this - * function simply passes through to - * the one implemented in the base class. - */ - virtual void fill_fe_values (const DoFHandler::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 &boundary) const; - - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_ansatz_points (const DoFHandler::cell_iterator &cell, - const Boundary &boundary, - vector > &ansatz_points) const; - - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_face_ansatz_points (const DoFHandler::face_iterator &face, - const Boundary &boundary, - vector > &ansatz_points) const; - - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_face_jacobians (const DoFHandler::face_iterator &face, - const Boundary &boundary, - const vector > &unit_points, - vector &face_jacobi_determinants) const; - - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_subface_jacobians (const DoFHandler::face_iterator &face, - const unsigned int subface_no, - const vector > &unit_points, - vector &face_jacobi_determinants) const; - - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_normal_vectors (const DoFHandler::cell_iterator &cell, - const unsigned int face_no, - const Boundary &boundary, - const vector > &unit_points, - vector > &normal_vectors) const; + Point linear_shape_grad(const unsigned int i, + const Point& p) const; +}; - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_normal_vectors (const DoFHandler::cell_iterator &cell, - const unsigned int subface_no, - const unsigned int face_no, - const vector > &unit_points, - vector > &normal_vectors) const; - /** - * Refer to the base class for detailed - * information on this function. - */ - virtual void get_local_mass_matrix (const DoFHandler::cell_iterator &cell, - const Boundary &boundary, - dFMatrix &local_mass_matrix) const; -}; diff --git a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc index 49ec2a0745..3fb43efa0c 100644 --- a/deal.II/deal.II/source/fe/fe_lib.quadratic.cc +++ b/deal.II/deal.II/source/fe/fe_lib.quadratic.cc @@ -15,7 +15,50 @@ template <> FEQuadraticSub<1>::FEQuadraticSub () : - FiniteElement<1> (1, 1) {}; + FiniteElement<1> (1, 1) { +/* + Get the prolongation matrices by the following little maple script: + + phi[0] := proc(xi) (1-xi)*(1-2*xi); end; + phi[1] := proc(xi) xi*(2*xi-1); end; + phi[2] := proc(xi) 4*xi*(1-xi); end; + + points[0] := array(0..2, [0, 1/2, 1/4]); + points[1] := array(0..2, [1/2, 1, 3/4]); + + prolongation := array(0..1,0..2, 0..2); + + for i from 0 to 1 do + for j from 0 to 2 do + for k from 0 to 2 do + prolongation[i,j,k] := phi[k](points[i][j]); + od; + od; + od; + + readlib(C); + C(prolongation); +*/ + + prolongation[0](0,0) = 1.0; + prolongation[0](0,1) = 0.0; + prolongation[0](0,2) = 0.0; + prolongation[0](1,0) = 0.0; + prolongation[0](1,1) = 0.0; + prolongation[0](1,2) = 1.0; + prolongation[0](2,0) = 3.0/8.0; + prolongation[0](2,1) = -1.0/8.0; + prolongation[0](2,2) = 3.0/4.0; + prolongation[1](0,0) = 0.0; + prolongation[1](0,1) = 0.0; + prolongation[1](0,2) = 1.0; + prolongation[1](1,0) = 0.0; + prolongation[1](1,1) = 1.0; + prolongation[1](1,2) = 0.0; + prolongation[1](2,0) = -1.0/8.0; + prolongation[1](2,1) = 3.0/8.0; + prolongation[1](2,2) = 3.0/4.0; +}; @@ -52,12 +95,29 @@ FEQuadraticSub<1>::shape_value(const unsigned int i, case 2: return 4*xi*(1-xi); } return 0.; -} +}; template <> inline +double +FEQuadraticSub<1>::linear_shape_value(const unsigned int i, + const Point<1> &p) const +{ + Assert((i<2), ExcInvalidIndex(i)); + const double xi = p(0); + switch (i) + { + case 0: return 1.-xi; + case 1: return xi; + } + return 0.; +}; + + + +template <> Point<1> FEQuadraticSub<1>::shape_grad(const unsigned int i, const Point<1> &p) const @@ -75,6 +135,23 @@ FEQuadraticSub<1>::shape_grad(const unsigned int i, +template <> +inline +Point<1> +FELinear<1>::linear_shape_grad(const unsigned int i, + const Point<1>&) const +{ + Assert((i<2), ExcInvalidIndex(i)); + switch (i) + { + case 0: return Point<1>(-1.); + case 1: return Point<1>(1.); + } + return Point<1>(); +}; + + + template <> void FEQuadraticSub<1>::get_ansatz_points (const typename DoFHandler<1>::cell_iterator &cell, const Boundary<1> &boundary, @@ -172,9 +249,369 @@ FEQuadraticSub<2>::FEQuadraticSub () : interface_constraints(2,1) = 3./8.; interface_constraints(2,2) = 3./4.; - // still implement restriction - // and prolongation - Assert (false, ExcNotImplemented()); +/* + Get the prolongation matrices by the following little maple script: + + phi[0] := proc(xi,eta) (1-xi)*( 2*xi-1) * (1-eta)*( 2*eta-1); end; + phi[1] := proc(xi,eta) xi *(-2*xi+1) * (1-eta)*( 2*eta-1); end; + phi[2] := proc(xi,eta) xi *(-2*xi+1) * eta *(-2*eta+1); end; + phi[3] := proc(xi,eta) (1-xi)*( 2*xi-1) * eta *(-2*eta+1); end; + phi[4] := proc(xi,eta) 4 * (1-xi)*xi * (1-eta)*(1-2*eta); end; + phi[5] := proc(xi,eta) 4 * xi *(-1+2*xi) * (1-eta)*eta; end; + phi[6] := proc(xi,eta) 4 * (1-xi)*xi * eta *(-1+2*eta);end; + phi[7] := proc(xi,eta) 4 * (1-xi)*(1-2*xi) * (1-eta)*eta; end; + phi[8] := proc(xi,eta) 16 * xi*(1-xi) * eta*(1-eta); end; + + points_x[0] := array(0..8, [0, 1/2, 1/2, 0, 1/4, 1/2, 1/4, 0, 1/4]); + points_y[0] := array(0..8, [0, 0, 1/2, 1/2, 0, 1/4, 1/2, 1/4, 1/4]); + + points_x[1] := array(0..8, [1/2, 1, 1, 1/2, 3/4, 1, 3/4, 1/2, 3/4]); + points_y[1] := array(0..8, [0, 0, 1/2, 1/2, 0, 1/4, 1/2, 1/4, 1/4]); + + points_x[2] := array(0..8, [1/2, 1, 1, 1/2, 3/4, 1, 3/4, 1/2, 3/4]); + points_y[2] := array(0..8, [1/2, 1/2, 1, 1, 1/2, 3/4, 1, 3/4, 3/4]); + + points_x[3] := array(0..8, [0, 1/2, 1/2, 0, 1/4, 1/2, 1/4, 0, 1/4]); + points_y[3] := array(0..8, [1/2, 1/2, 1, 1, 1/2, 3/4, 1, 3/4, 3/4]); + + prolongation := array(0..3,0..8, 0..8); + + for i from 0 to 3 do + for j from 0 to 8 do + for k from 0 to 8 do + prolongation[i,j,k] := phi[k](points_x[i][j], points_y[i][j]); + od; + od; + od; + + readlib(C); + C(prolongation); +*/ + + prolongation[0](0,0) = 1.0; + prolongation[0](0,1) = 0.0; + prolongation[0](0,2) = 0.0; + prolongation[0](0,3) = 0.0; + prolongation[0](0,4) = 0.0; + prolongation[0](0,5) = 0.0; + prolongation[0](0,6) = 0.0; + prolongation[0](0,7) = 0.0; + prolongation[0](0,8) = 0.0; + prolongation[0](1,0) = 0.0; + prolongation[0](1,1) = 0.0; + prolongation[0](1,2) = 0.0; + prolongation[0](1,3) = 0.0; + prolongation[0](1,4) = 1.0; + prolongation[0](1,5) = 0.0; + prolongation[0](1,6) = 0.0; + prolongation[0](1,7) = 0.0; + prolongation[0](1,8) = 0.0; + prolongation[0](2,0) = 0.0; + prolongation[0](2,1) = 0.0; + prolongation[0](2,2) = 0.0; + prolongation[0](2,3) = 0.0; + prolongation[0](2,4) = 0.0; + prolongation[0](2,5) = 0.0; + prolongation[0](2,6) = 0.0; + prolongation[0](2,7) = 0.0; + prolongation[0](2,8) = 1.0; + prolongation[0](3,0) = 0.0; + prolongation[0](3,1) = 0.0; + prolongation[0](3,2) = 0.0; + prolongation[0](3,3) = 0.0; + prolongation[0](3,4) = 0.0; + prolongation[0](3,5) = 0.0; + prolongation[0](3,6) = 0.0; + prolongation[0](3,7) = 1.0; + prolongation[0](3,8) = 0.0; + prolongation[0](4,0) = 3.0/8.0; + prolongation[0](4,1) = -1.0/8.0; + prolongation[0](4,2) = 0.0; + prolongation[0](4,3) = 0.0; + prolongation[0](4,4) = 3.0/4.0; + prolongation[0](4,5) = 0.0; + prolongation[0](4,6) = 0.0; + prolongation[0](4,7) = 0.0; + prolongation[0](4,8) = 0.0; + prolongation[0](5,0) = 0.0; + prolongation[0](5,1) = 0.0; + prolongation[0](5,2) = 0.0; + prolongation[0](5,3) = 0.0; + prolongation[0](5,4) = 3.0/8.0; + prolongation[0](5,5) = 0.0; + prolongation[0](5,6) = -1.0/8.0; + prolongation[0](5,7) = 0.0; + prolongation[0](5,8) = 3.0/4.0; + prolongation[0](6,0) = 0.0; + prolongation[0](6,1) = 0.0; + prolongation[0](6,2) = 0.0; + prolongation[0](6,3) = 0.0; + prolongation[0](6,4) = 0.0; + prolongation[0](6,5) = -1.0/8.0; + prolongation[0](6,6) = 0.0; + prolongation[0](6,7) = 3.0/8.0; + prolongation[0](6,8) = 3.0/4.0; + prolongation[0](7,0) = 3.0/8.0; + prolongation[0](7,1) = 0.0; + prolongation[0](7,2) = 0.0; + prolongation[0](7,3) = -1.0/8.0; + prolongation[0](7,4) = 0.0; + prolongation[0](7,5) = 0.0; + prolongation[0](7,6) = 0.0; + prolongation[0](7,7) = 3.0/4.0; + prolongation[0](7,8) = 0.0; + prolongation[0](8,0) = 9.0/64.0; + prolongation[0](8,1) = -3.0/64.0; + prolongation[0](8,2) = 1.0/64.0; + prolongation[0](8,3) = -3.0/64.0; + prolongation[0](8,4) = 9.0/32.0; + prolongation[0](8,5) = -3.0/32.0; + prolongation[0](8,6) = -3.0/32.0; + prolongation[0](8,7) = 9.0/32.0; + prolongation[0](8,8) = 9.0/16.0; + prolongation[1](0,0) = 0.0; + prolongation[1](0,1) = 0.0; + prolongation[1](0,2) = 0.0; + prolongation[1](0,3) = 0.0; + prolongation[1](0,4) = 1.0; + prolongation[1](0,5) = 0.0; + prolongation[1](0,6) = 0.0; + prolongation[1](0,7) = 0.0; + prolongation[1](0,8) = 0.0; + prolongation[1](1,0) = 0.0; + prolongation[1](1,1) = 1.0; + prolongation[1](1,2) = 0.0; + prolongation[1](1,3) = 0.0; + prolongation[1](1,4) = 0.0; + prolongation[1](1,5) = 0.0; + prolongation[1](1,6) = 0.0; + prolongation[1](1,7) = 0.0; + prolongation[1](1,8) = 0.0; + prolongation[1](2,0) = 0.0; + prolongation[1](2,1) = 0.0; + prolongation[1](2,2) = 0.0; + prolongation[1](2,3) = 0.0; + prolongation[1](2,4) = 0.0; + prolongation[1](2,5) = 1.0; + prolongation[1](2,6) = 0.0; + prolongation[1](2,7) = 0.0; + prolongation[1](2,8) = 0.0; + prolongation[1](3,0) = 0.0; + prolongation[1](3,1) = 0.0; + prolongation[1](3,2) = 0.0; + prolongation[1](3,3) = 0.0; + prolongation[1](3,4) = 0.0; + prolongation[1](3,5) = 0.0; + prolongation[1](3,6) = 0.0; + prolongation[1](3,7) = 0.0; + prolongation[1](3,8) = 1.0; + prolongation[1](4,0) = -1.0/8.0; + prolongation[1](4,1) = 3.0/8.0; + prolongation[1](4,2) = 0.0; + prolongation[1](4,3) = 0.0; + prolongation[1](4,4) = 3.0/4.0; + prolongation[1](4,5) = 0.0; + prolongation[1](4,6) = 0.0; + prolongation[1](4,7) = 0.0; + prolongation[1](4,8) = 0.0; + prolongation[1](5,0) = 0.0; + prolongation[1](5,1) = 3.0/8.0; + prolongation[1](5,2) = -1.0/8.0; + prolongation[1](5,3) = 0.0; + prolongation[1](5,4) = 0.0; + prolongation[1](5,5) = 3.0/4.0; + prolongation[1](5,6) = 0.0; + prolongation[1](5,7) = 0.0; + prolongation[1](5,8) = 0.0; + prolongation[1](6,0) = 0.0; + prolongation[1](6,1) = 0.0; + prolongation[1](6,2) = 0.0; + prolongation[1](6,3) = 0.0; + prolongation[1](6,4) = 0.0; + prolongation[1](6,5) = 3.0/8.0; + prolongation[1](6,6) = 0.0; + prolongation[1](6,7) = -1.0/8.0; + prolongation[1](6,8) = 3.0/4.0; + prolongation[1](7,0) = 0.0; + prolongation[1](7,1) = 0.0; + prolongation[1](7,2) = 0.0; + prolongation[1](7,3) = 0.0; + prolongation[1](7,4) = 3.0/8.0; + prolongation[1](7,5) = 0.0; + prolongation[1](7,6) = -1.0/8.0; + prolongation[1](7,7) = 0.0; + prolongation[1](7,8) = 3.0/4.0; + prolongation[1](8,0) = -3.0/64.0; + prolongation[1](8,1) = 9.0/64.0; + prolongation[1](8,2) = -3.0/64.0; + prolongation[1](8,3) = 1.0/64.0; + prolongation[1](8,4) = 9.0/32.0; + prolongation[1](8,5) = 9.0/32.0; + prolongation[1](8,6) = -3.0/32.0; + prolongation[1](8,7) = -3.0/32.0; + prolongation[1](8,8) = 9.0/16.0; + prolongation[2](0,0) = 0.0; + prolongation[2](0,1) = 0.0; + prolongation[2](0,2) = 0.0; + prolongation[2](0,3) = 0.0; + prolongation[2](0,4) = 0.0; + prolongation[2](0,5) = 0.0; + prolongation[2](0,6) = 0.0; + prolongation[2](0,7) = 0.0; + prolongation[2](0,8) = 1.0; + prolongation[2](1,0) = 0.0; + prolongation[2](1,1) = 0.0; + prolongation[2](1,2) = 0.0; + prolongation[2](1,3) = 0.0; + prolongation[2](1,4) = 0.0; + prolongation[2](1,5) = 1.0; + prolongation[2](1,6) = 0.0; + prolongation[2](1,7) = 0.0; + prolongation[2](1,8) = 0.0; + prolongation[2](2,0) = 0.0; + prolongation[2](2,1) = 0.0; + prolongation[2](2,2) = 1.0; + prolongation[2](2,3) = 0.0; + prolongation[2](2,4) = 0.0; + prolongation[2](2,5) = 0.0; + prolongation[2](2,6) = 0.0; + prolongation[2](2,7) = 0.0; + prolongation[2](2,8) = 0.0; + prolongation[2](3,0) = 0.0; + prolongation[2](3,1) = 0.0; + prolongation[2](3,2) = 0.0; + prolongation[2](3,3) = 0.0; + prolongation[2](3,4) = 0.0; + prolongation[2](3,5) = 0.0; + prolongation[2](3,6) = 1.0; + prolongation[2](3,7) = 0.0; + prolongation[2](3,8) = 0.0; + prolongation[2](4,0) = 0.0; + prolongation[2](4,1) = 0.0; + prolongation[2](4,2) = 0.0; + prolongation[2](4,3) = 0.0; + prolongation[2](4,4) = 0.0; + prolongation[2](4,5) = 3.0/8.0; + prolongation[2](4,6) = 0.0; + prolongation[2](4,7) = -1.0/8.0; + prolongation[2](4,8) = 3.0/4.0; + prolongation[2](5,0) = 0.0; + prolongation[2](5,1) = -1.0/8.0; + prolongation[2](5,2) = 3.0/8.0; + prolongation[2](5,3) = 0.0; + prolongation[2](5,4) = 0.0; + prolongation[2](5,5) = 3.0/4.0; + prolongation[2](5,6) = 0.0; + prolongation[2](5,7) = 0.0; + prolongation[2](5,8) = 0.0; + prolongation[2](6,0) = 0.0; + prolongation[2](6,1) = 0.0; + prolongation[2](6,2) = 3.0/8.0; + prolongation[2](6,3) = -1.0/8.0; + prolongation[2](6,4) = 0.0; + prolongation[2](6,5) = 0.0; + prolongation[2](6,6) = 3.0/4.0; + prolongation[2](6,7) = 0.0; + prolongation[2](6,8) = 0.0; + prolongation[2](7,0) = 0.0; + prolongation[2](7,1) = 0.0; + prolongation[2](7,2) = 0.0; + prolongation[2](7,3) = 0.0; + prolongation[2](7,4) = -1.0/8.0; + prolongation[2](7,5) = 0.0; + prolongation[2](7,6) = 3.0/8.0; + prolongation[2](7,7) = 0.0; + prolongation[2](7,8) = 3.0/4.0; + prolongation[2](8,0) = 1.0/64.0; + prolongation[2](8,1) = -3.0/64.0; + prolongation[2](8,2) = 9.0/64.0; + prolongation[2](8,3) = -3.0/64.0; + prolongation[2](8,4) = -3.0/32.0; + prolongation[2](8,5) = 9.0/32.0; + prolongation[2](8,6) = 9.0/32.0; + prolongation[2](8,7) = -3.0/32.0; + prolongation[2](8,8) = 9.0/16.0; + prolongation[3](0,0) = 0.0; + prolongation[3](0,1) = 0.0; + prolongation[3](0,2) = 0.0; + prolongation[3](0,3) = 0.0; + prolongation[3](0,4) = 0.0; + prolongation[3](0,5) = 0.0; + prolongation[3](0,6) = 0.0; + prolongation[3](0,7) = 1.0; + prolongation[3](0,8) = 0.0; + prolongation[3](1,0) = 0.0; + prolongation[3](1,1) = 0.0; + prolongation[3](1,2) = 0.0; + prolongation[3](1,3) = 0.0; + prolongation[3](1,4) = 0.0; + prolongation[3](1,5) = 0.0; + prolongation[3](1,6) = 0.0; + prolongation[3](1,7) = 0.0; + prolongation[3](1,8) = 1.0; + prolongation[3](2,0) = 0.0; + prolongation[3](2,1) = 0.0; + prolongation[3](2,2) = 0.0; + prolongation[3](2,3) = 0.0; + prolongation[3](2,4) = 0.0; + prolongation[3](2,5) = 0.0; + prolongation[3](2,6) = 1.0; + prolongation[3](2,7) = 0.0; + prolongation[3](2,8) = 0.0; + prolongation[3](3,0) = 0.0; + prolongation[3](3,1) = 0.0; + prolongation[3](3,2) = 0.0; + prolongation[3](3,3) = 1.0; + prolongation[3](3,4) = 0.0; + prolongation[3](3,5) = 0.0; + prolongation[3](3,6) = 0.0; + prolongation[3](3,7) = 0.0; + prolongation[3](3,8) = 0.0; + prolongation[3](4,0) = 0.0; + prolongation[3](4,1) = 0.0; + prolongation[3](4,2) = 0.0; + prolongation[3](4,3) = 0.0; + prolongation[3](4,4) = 0.0; + prolongation[3](4,5) = -1.0/8.0; + prolongation[3](4,6) = 0.0; + prolongation[3](4,7) = 3.0/8.0; + prolongation[3](4,8) = 3.0/4.0; + prolongation[3](5,0) = 0.0; + prolongation[3](5,1) = 0.0; + prolongation[3](5,2) = 0.0; + prolongation[3](5,3) = 0.0; + prolongation[3](5,4) = -1.0/8.0; + prolongation[3](5,5) = 0.0; + prolongation[3](5,6) = 3.0/8.0; + prolongation[3](5,7) = 0.0; + prolongation[3](5,8) = 3.0/4.0; + prolongation[3](6,0) = 0.0; + prolongation[3](6,1) = 0.0; + prolongation[3](6,2) = -1.0/8.0; + prolongation[3](6,3) = 3.0/8.0; + prolongation[3](6,4) = 0.0; + prolongation[3](6,5) = 0.0; + prolongation[3](6,6) = 3.0/4.0; + prolongation[3](6,7) = 0.0; + prolongation[3](6,8) = 0.0; + prolongation[3](7,0) = -1.0/8.0; + prolongation[3](7,1) = 0.0; + prolongation[3](7,2) = 0.0; + prolongation[3](7,3) = 3.0/8.0; + prolongation[3](7,4) = 0.0; + prolongation[3](7,5) = 0.0; + prolongation[3](7,6) = 0.0; + prolongation[3](7,7) = 3.0/4.0; + prolongation[3](7,8) = 0.0; + prolongation[3](8,0) = -3.0/64.0; + prolongation[3](8,1) = 1.0/64.0; + prolongation[3](8,2) = -3.0/64.0; + prolongation[3](8,3) = 9.0/64.0; + prolongation[3](8,4) = -3.0/32.0; + prolongation[3](8,5) = -3.0/32.0; + prolongation[3](8,6) = 9.0/32.0; + prolongation[3](8,7) = 9.0/32.0; + prolongation[3](8,8) = 9.0/16.0; }; @@ -678,15 +1115,15 @@ void FEQuadraticSub<2>::get_normal_vectors (const DoFHandler<2>::cell_iterator & 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 { +void FEQuadraticSub::fill_fe_values (const DoFHandler::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 &boundary) const { Assert (jacobians.size() == unit_points.size(), ExcWrongFieldDimension(jacobians.size(), unit_points.size())); Assert (q_points.size() == unit_points.size(), @@ -694,230 +1131,79 @@ void FEQuadraticSub::fill_fe_values (const DoFHandler::cell_iterator & 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 ()); + + unsigned int n_points=unit_points.size(); - Assert (false, ExcNotImplemented()); -}; + Point vertices[GeometryInfo::vertices_per_cell]; + for (unsigned int l=0; l::vertices_per_cell; ++l) + vertices[l] = cell->vertex(l); + + if (compute_q_points) + { + // initialize points to zero + for (unsigned int i=0; i (); + + // note: let x_l be the vector of the + // lth quadrature point in real space and + // xi_l that on the unit cell, let further + // p_j be the vector of the jth vertex + // of the cell in real space and + // N_j(xi_l) be the value of the associated + // basis function at xi_l, then + // x_l(xi_l) = sum_j p_j N_j(xi_l) + // + // Here, N_j is the *linear* basis function, + // not that of the finite element, since we + // use a subparametric mapping + for (unsigned int j=0; j::vertices_per_cell; ++j) + for (unsigned int l=0; l::vertices_per_cell; ++s) + { + // we want the linear transform, + // so use that function + const Point gradient = linear_shape_grad (s, unit_points[l]); + for (unsigned int i=0; i -void FECubic::get_local_mass_matrix (const DoFHandler::cell_iterator &, - const Boundary &, - dFMatrix &) const { - Assert (false, ExcNotImplemented()); + // compute ansatz points, which are + // the corners for linear elements + if (compute_ansatz_points) + get_ansatz_points (cell, boundary, ansatz_points); }; -- 2.39.5