From d9a950ba736ca8210dabbc33d6988eac0138cefc Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 23 Apr 2003 22:59:17 +0000 Subject: [PATCH] Single out the functions that initialize matrices from the constructor, to make the whole process of initializing finite element classes a little more transparent. Also makes it simple to replace precomputed matrices with something that computes these values on the fly. git-svn-id: https://svn.dealii.org/trunk@7463 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_q.h | 27 +- deal.II/deal.II/source/fe/fe_q.cc | 1646 +++++++++++++++-------------- 2 files changed, 865 insertions(+), 808 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index cb4ec781fd..218949d220 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -624,6 +624,9 @@ class FE_Q : public FiniteElement * chosen to weigh the simplicity * aspect a little more than * proper design. + * + * This function is called from + * the constructor. */ static std::vector @@ -633,12 +636,34 @@ class FE_Q : public FiniteElement /** * This is an analogon to the * previous function, but working - * on faces. + * on faces. Called from the + * constructor. */ static std::vector face_lexicographic_to_hierarchic_numbering (const unsigned int degree); + /** + * Initialize the hanging node + * constraints matrices. Called + * from the constructor. + */ + void initialize_constraints (); + + /** + * Initialize the embedding + * matrices. Called from the + * constructor. + */ + void initialize_embedding (); + + /** + * Initialize the restriction + * matrices. Called from the + * constructor. + */ + void initialize_restriction (); + /** * Initialize the * @p{unit_support_points} field diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index b39b401616..8a85752502 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -51,779 +51,213 @@ FE_Q::FE_Q (const unsigned int degree) face_renumber(face_lexicographic_to_hierarchic_numbering (degree)), polynomial_space(Polynomials::LagrangeEquidistant::generate_complete_basis(degree)) { - // copy constraint matrices if they - // are defined. otherwise leave them - // at invalid size - if ((dim>1) && (degree +FiniteElement * +FE_Q::clone() const +{ + return new FE_Q(degree); +} + + + +template +double +FE_Q::shape_value (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + return polynomial_space.compute_value(renumber_inverse[i], p); +} + + +template +double +FE_Q::shape_value_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + Assert (component == 0, ExcIndexRange (component, 0, 1)); + return polynomial_space.compute_value(renumber_inverse[i], p); +} + + + +template +Tensor<1,dim> +FE_Q::shape_grad (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + return polynomial_space.compute_grad(renumber_inverse[i], p); +} + + + +template +Tensor<1,dim> +FE_Q::shape_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + Assert (component == 0, ExcIndexRange (component, 0, 1)); + return polynomial_space.compute_grad(renumber_inverse[i], p); +} + + + +template +Tensor<2,dim> +FE_Q::shape_grad_grad (const unsigned int i, + const Point &p) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + return polynomial_space.compute_grad_grad(renumber_inverse[i], p); +} + + + +template +Tensor<2,dim> +FE_Q::shape_grad_grad_component (const unsigned int i, + const Point &p, + const unsigned int component) const +{ + Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); + Assert (component == 0, ExcIndexRange (component, 0, 1)); + return polynomial_space.compute_grad_grad(renumber_inverse[i], p); +} + + +//---------------------------------------------------------------------- +// Auxiliary functions +//---------------------------------------------------------------------- + + + +template +void FE_Q::initialize_unit_support_points () +{ + // number of points: (degree+1)^dim + unsigned int n = degree+1; + for (unsigned int i=1; iunit_support_points.resize(n); + + const double step = 1./degree; + Point p; + + unsigned int k=0; + for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz) + for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy) + for (unsigned int ix=0; ix<=degree; ++ix) + { + p(0) = ix * step; + if (dim>1) + p(1) = iy * step; + if (dim>2) + p(2) = iz * step; + + this->unit_support_points[renumber[k++]] = p; + }; +} + + +#if deal_II_dimension == 1 + +template <> +void FE_Q<1>::initialize_unit_face_support_points () +{ + // no faces in 1d, so nothing to do +} + +#endif + + +template +void FE_Q::initialize_unit_face_support_points () +{ + const unsigned int codim = dim-1; + + // number of points: (degree+1)^codim + unsigned int n = degree+1; + for (unsigned int i=1; iunit_face_support_points.resize(n); + + const double step = 1./degree; + Point p; + + unsigned int k=0; + for (unsigned int iz=0; iz <= ((codim>2) ? degree : 0) ; ++iz) + for (unsigned int iy=0; iy <= ((codim>1) ? degree : 0) ; ++iy) + for (unsigned int ix=0; ix<=degree; ++ix) + { + p(0) = ix * step; + if (codim>1) + p(1) = iy * step; + if (codim>2) + p(2) = iz * step; + + this->unit_face_support_points[face_renumber[k++]] = p; + }; +} + + + +template +std::vector +FE_Q::get_dpo_vector(const unsigned int deg) +{ + std::vector dpo(dim+1, static_cast(1)); + for (unsigned int i=1; i +std::vector +FE_Q::lexicographic_to_hierarchic_numbering (const FiniteElementData &fe_data, + const unsigned int degree) +{ + std::vector renumber (fe_data.dofs_per_cell); + + const unsigned int n = degree+1; + + + if (degree == 0) { - this->interface_constraints. - TableBase<2,double>::reinit (this->interface_constraints_size()); - this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]); + Assert ((fe_data.dofs_per_vertex == 0) && + ((fe_data.dofs_per_line == 0) || (dim == 1)) && + ((fe_data.dofs_per_quad == 0) || (dim == 2)) && + ((fe_data.dofs_per_hex == 0) || (dim == 3)), + ExcInternalError()); + renumber[0] = 0; }; - // next copy over embedding - // matrices if they are defined - if ((degree < Matrices::n_embedding_matrices+1) && - (Matrices::embedding[degree-1][0] != 0)) - for (unsigned int c=0; c::children_per_cell; ++c) - { - this->prolongation[c].reinit (this->dofs_per_cell, - this->dofs_per_cell); - this->prolongation[c].fill (Matrices::embedding[degree-1][c]); - }; - - // then fill restriction - // matrices. they are hardcoded for - // the first few elements. in - // contrast to the other matrices, - // these are not stored in the - // files fe_q_[123]d.cc, since they - // contain only a rather small - // number of zeros, and storing - // them element-wise is more - // expensive than just setting the - // nonzero elements as done here - for (unsigned int c=0; c::children_per_cell; ++c) - this->restriction[c].reinit (this->dofs_per_cell, this->dofs_per_cell); - switch (dim) - { - case 1: // 1d - { - switch (degree) - { - case 1: - this->restriction[0](0,0) = 1; - this->restriction[1](1,1) = 1; - break; - case 2: - this->restriction[0](0,0) = 1; - this->restriction[0](2,1) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](1,1) = 1; - break; - case 3: - this->restriction[0](0,0) = 1; - this->restriction[0](2,3) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](3,2) = 1; - break; - case 4: - this->restriction[0](0,0) = 1; - this->restriction[0](2,3) = 1; - this->restriction[0](3,1) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](3,0) = 1; - this->restriction[1](4,3) = 1; - break; - - default: - { - // in case we don't - // have the matrices - // (yet), reset them to - // zero size. this does - // not prevent the use - // of this FE, but will - // prevent the use of - // these matrices - for (unsigned int i=0; - i::children_per_cell; - ++i) - this->restriction[i].reinit(0,0); - }; - } - break; - }; - - case 2: // 2d - { - switch (degree) - { - case 1: - this->restriction[0](0,0) = 1; - this->restriction[1](1,1) = 1; - this->restriction[2](2,2) = 1; - this->restriction[3](3,3) = 1; - break; - case 2: - this->restriction[0](0,0) = 1; - this->restriction[0](4,1) = 1; - this->restriction[0](7,3) = 1; - this->restriction[0](8,2) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](4,0) = 1; - this->restriction[1](5,2) = 1; - this->restriction[1](8,3) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](5,1) = 1; - this->restriction[2](6,3) = 1; - this->restriction[2](8,0) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](6,2) = 1; - this->restriction[3](7,0) = 1; - this->restriction[3](8,1) = 1; - break; - case 3: - this->restriction[0](0,0) = 1; - this->restriction[0](4,5) = 1; - this->restriction[0](10,11) = 1; - this->restriction[0](12,15) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](5,4) = 1; - this->restriction[1](6,7) = 1; - this->restriction[1](13,14) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](7,6) = 1; - this->restriction[2](9,8) = 1; - this->restriction[2](15,12) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](8,9) = 1; - this->restriction[3](11,10) = 1; - this->restriction[3](14,13) = 1; - break; - case 4: - this->restriction[0](0,0) = 1; - this->restriction[0](4,5) = 1; - this->restriction[0](5,1) = 1; - this->restriction[0](13,14) = 1; - this->restriction[0](14,3) = 1; - this->restriction[0](16,20) = 1; - this->restriction[0](17,8) = 1; - this->restriction[0](19,11) = 1; - this->restriction[0](20,2) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](5,0) = 1; - this->restriction[1](6,5) = 1; - this->restriction[1](7,8) = 1; - this->restriction[1](8,2) = 1; - this->restriction[1](17,14) = 1; - this->restriction[1](18,20) = 1; - this->restriction[1](20,3) = 1; - this->restriction[1](21,11) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](8,1) = 1; - this->restriction[2](9,8) = 1; - this->restriction[2](11,3) = 1; - this->restriction[2](12,11) = 1; - this->restriction[2](20,0) = 1; - this->restriction[2](21,5) = 1; - this->restriction[2](23,14) = 1; - this->restriction[2](24,20) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](10,11) = 1; - this->restriction[3](11,2) = 1; - this->restriction[3](14,0) = 1; - this->restriction[3](15,14) = 1; - this->restriction[3](19,5) = 1; - this->restriction[3](20,1) = 1; - this->restriction[3](22,20) = 1; - this->restriction[3](23,8) = 1; - break; - - default: - { - // in case we don't - // have the matrices - // (yet), reset them to - // zero size. this does - // not prevent the use - // of this FE, but will - // prevent the use of - // these matrices - for (unsigned int i=0; - i::children_per_cell; - ++i) - this->restriction[i].reinit(0,0); - }; - } - break; - }; - - case 3: // 3d - { - switch (degree) - { - case 1: - this->restriction[0](0,0) = 1; - this->restriction[1](1,1) = 1; - this->restriction[2](2,2) = 1; - this->restriction[3](3,3) = 1; - this->restriction[4](4,4) = 1; - this->restriction[5](5,5) = 1; - this->restriction[6](6,6) = 1; - this->restriction[7](7,7) = 1; - break; - case 2: - this->restriction[0](0,0) = 1; - this->restriction[0](8,1) = 1; - this->restriction[0](11,3) = 1; - this->restriction[0](16,4) = 1; - this->restriction[0](20,2) = 1; - this->restriction[0](22,5) = 1; - this->restriction[0](25,7) = 1; - this->restriction[0](26,6) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](8,0) = 1; - this->restriction[1](9,2) = 1; - this->restriction[1](17,5) = 1; - this->restriction[1](20,3) = 1; - this->restriction[1](22,4) = 1; - this->restriction[1](23,6) = 1; - this->restriction[1](26,7) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](9,1) = 1; - this->restriction[2](10,3) = 1; - this->restriction[2](18,6) = 1; - this->restriction[2](20,0) = 1; - this->restriction[2](23,5) = 1; - this->restriction[2](24,7) = 1; - this->restriction[2](26,4) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](10,2) = 1; - this->restriction[3](11,0) = 1; - this->restriction[3](19,7) = 1; - this->restriction[3](20,1) = 1; - this->restriction[3](24,6) = 1; - this->restriction[3](25,4) = 1; - this->restriction[3](26,5) = 1; - this->restriction[4](4,4) = 1; - this->restriction[4](12,5) = 1; - this->restriction[4](15,7) = 1; - this->restriction[4](16,0) = 1; - this->restriction[4](21,6) = 1; - this->restriction[4](22,1) = 1; - this->restriction[4](25,3) = 1; - this->restriction[4](26,2) = 1; - this->restriction[5](5,5) = 1; - this->restriction[5](12,4) = 1; - this->restriction[5](13,6) = 1; - this->restriction[5](17,1) = 1; - this->restriction[5](21,7) = 1; - this->restriction[5](22,0) = 1; - this->restriction[5](23,2) = 1; - this->restriction[5](26,3) = 1; - this->restriction[6](6,6) = 1; - this->restriction[6](13,5) = 1; - this->restriction[6](14,7) = 1; - this->restriction[6](18,2) = 1; - this->restriction[6](21,4) = 1; - this->restriction[6](23,1) = 1; - this->restriction[6](24,3) = 1; - this->restriction[6](26,0) = 1; - this->restriction[7](7,7) = 1; - this->restriction[7](14,6) = 1; - this->restriction[7](15,4) = 1; - this->restriction[7](19,3) = 1; - this->restriction[7](21,5) = 1; - this->restriction[7](24,2) = 1; - this->restriction[7](25,0) = 1; - this->restriction[7](26,1) = 1; - break; - case 3: - this->restriction[0](0,0) = 1; - this->restriction[0](8,9) = 1; - this->restriction[0](14,15) = 1; - this->restriction[0](24,25) = 1; - this->restriction[0](32,35) = 1; - this->restriction[0](40,43) = 1; - this->restriction[0](52,55) = 1; - this->restriction[0](56,63) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](9,8) = 1; - this->restriction[1](10,11) = 1; - this->restriction[1](26,27) = 1; - this->restriction[1](33,34) = 1; - this->restriction[1](41,42) = 1; - this->restriction[1](44,47) = 1; - this->restriction[1](57,62) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](11,10) = 1; - this->restriction[2](13,12) = 1; - this->restriction[2](28,29) = 1; - this->restriction[2](35,32) = 1; - this->restriction[2](46,45) = 1; - this->restriction[2](49,50) = 1; - this->restriction[2](61,58) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](12,13) = 1; - this->restriction[3](15,14) = 1; - this->restriction[3](30,31) = 1; - this->restriction[3](34,33) = 1; - this->restriction[3](48,51) = 1; - this->restriction[3](54,53) = 1; - this->restriction[3](60,59) = 1; - this->restriction[4](4,4) = 1; - this->restriction[4](16,17) = 1; - this->restriction[4](22,23) = 1; - this->restriction[4](25,24) = 1; - this->restriction[4](36,39) = 1; - this->restriction[4](42,41) = 1; - this->restriction[4](53,54) = 1; - this->restriction[4](58,61) = 1; - this->restriction[5](5,5) = 1; - this->restriction[5](17,16) = 1; - this->restriction[5](18,19) = 1; - this->restriction[5](27,26) = 1; - this->restriction[5](37,38) = 1; - this->restriction[5](43,40) = 1; - this->restriction[5](45,46) = 1; - this->restriction[5](59,60) = 1; - this->restriction[6](6,6) = 1; - this->restriction[6](19,18) = 1; - this->restriction[6](21,20) = 1; - this->restriction[6](29,28) = 1; - this->restriction[6](39,36) = 1; - this->restriction[6](47,44) = 1; - this->restriction[6](51,48) = 1; - this->restriction[6](63,56) = 1; - this->restriction[7](7,7) = 1; - this->restriction[7](20,21) = 1; - this->restriction[7](23,22) = 1; - this->restriction[7](31,30) = 1; - this->restriction[7](38,37) = 1; - this->restriction[7](50,49) = 1; - this->restriction[7](55,52) = 1; - this->restriction[7](62,57) = 1; - break; - case 4: - this->restriction[0](0,0) = 1; - this->restriction[0](8,9) = 1; - this->restriction[0](9,1) = 1; - this->restriction[0](17,18) = 1; - this->restriction[0](18,3) = 1; - this->restriction[0](32,33) = 1; - this->restriction[0](33,4) = 1; - this->restriction[0](44,48) = 1; - this->restriction[0](45,12) = 1; - this->restriction[0](47,15) = 1; - this->restriction[0](48,2) = 1; - this->restriction[0](62,66) = 1; - this->restriction[0](63,36) = 1; - this->restriction[0](65,21) = 1; - this->restriction[0](66,5) = 1; - this->restriction[0](89,93) = 1; - this->restriction[0](90,30) = 1; - this->restriction[0](92,42) = 1; - this->restriction[0](93,7) = 1; - this->restriction[0](98,111) = 1; - this->restriction[0](99,75) = 1; - this->restriction[0](101,57) = 1; - this->restriction[0](102,24) = 1; - this->restriction[0](107,84) = 1; - this->restriction[0](108,39) = 1; - this->restriction[0](110,27) = 1; - this->restriction[0](111,6) = 1; - this->restriction[1](1,1) = 1; - this->restriction[1](9,0) = 1; - this->restriction[1](10,9) = 1; - this->restriction[1](11,12) = 1; - this->restriction[1](12,2) = 1; - this->restriction[1](35,36) = 1; - this->restriction[1](36,5) = 1; - this->restriction[1](45,18) = 1; - this->restriction[1](46,48) = 1; - this->restriction[1](48,3) = 1; - this->restriction[1](49,15) = 1; - this->restriction[1](63,33) = 1; - this->restriction[1](64,66) = 1; - this->restriction[1](66,4) = 1; - this->restriction[1](67,21) = 1; - this->restriction[1](71,75) = 1; - this->restriction[1](72,24) = 1; - this->restriction[1](74,39) = 1; - this->restriction[1](75,6) = 1; - this->restriction[1](99,93) = 1; - this->restriction[1](100,111) = 1; - this->restriction[1](102,30) = 1; - this->restriction[1](103,57) = 1; - this->restriction[1](108,42) = 1; - this->restriction[1](109,84) = 1; - this->restriction[1](111,7) = 1; - this->restriction[1](112,27) = 1; - this->restriction[2](2,2) = 1; - this->restriction[2](12,1) = 1; - this->restriction[2](13,12) = 1; - this->restriction[2](15,3) = 1; - this->restriction[2](16,15) = 1; - this->restriction[2](38,39) = 1; - this->restriction[2](39,6) = 1; - this->restriction[2](48,0) = 1; - this->restriction[2](49,9) = 1; - this->restriction[2](51,18) = 1; - this->restriction[2](52,48) = 1; - this->restriction[2](74,36) = 1; - this->restriction[2](75,5) = 1; - this->restriction[2](77,75) = 1; - this->restriction[2](78,24) = 1; - this->restriction[2](81,42) = 1; - this->restriction[2](82,84) = 1; - this->restriction[2](84,7) = 1; - this->restriction[2](85,27) = 1; - this->restriction[2](108,33) = 1; - this->restriction[2](109,66) = 1; - this->restriction[2](111,4) = 1; - this->restriction[2](112,21) = 1; - this->restriction[2](117,93) = 1; - this->restriction[2](118,111) = 1; - this->restriction[2](120,30) = 1; - this->restriction[2](121,57) = 1; - this->restriction[3](3,3) = 1; - this->restriction[3](14,15) = 1; - this->restriction[3](15,2) = 1; - this->restriction[3](18,0) = 1; - this->restriction[3](19,18) = 1; - this->restriction[3](41,42) = 1; - this->restriction[3](42,7) = 1; - this->restriction[3](47,9) = 1; - this->restriction[3](48,1) = 1; - this->restriction[3](50,48) = 1; - this->restriction[3](51,12) = 1; - this->restriction[3](80,84) = 1; - this->restriction[3](81,39) = 1; - this->restriction[3](83,27) = 1; - this->restriction[3](84,6) = 1; - this->restriction[3](92,33) = 1; - this->restriction[3](93,4) = 1; - this->restriction[3](95,93) = 1; - this->restriction[3](96,30) = 1; - this->restriction[3](107,66) = 1; - this->restriction[3](108,36) = 1; - this->restriction[3](110,21) = 1; - this->restriction[3](111,5) = 1; - this->restriction[3](116,111) = 1; - this->restriction[3](117,75) = 1; - this->restriction[3](119,57) = 1; - this->restriction[3](120,24) = 1; - this->restriction[4](4,4) = 1; - this->restriction[4](20,21) = 1; - this->restriction[4](21,5) = 1; - this->restriction[4](29,30) = 1; - this->restriction[4](30,7) = 1; - this->restriction[4](33,0) = 1; - this->restriction[4](34,33) = 1; - this->restriction[4](53,57) = 1; - this->restriction[4](54,24) = 1; - this->restriction[4](56,27) = 1; - this->restriction[4](57,6) = 1; - this->restriction[4](65,9) = 1; - this->restriction[4](66,1) = 1; - this->restriction[4](68,66) = 1; - this->restriction[4](69,36) = 1; - this->restriction[4](90,18) = 1; - this->restriction[4](91,93) = 1; - this->restriction[4](93,3) = 1; - this->restriction[4](94,42) = 1; - this->restriction[4](101,48) = 1; - this->restriction[4](102,12) = 1; - this->restriction[4](104,111) = 1; - this->restriction[4](105,75) = 1; - this->restriction[4](110,15) = 1; - this->restriction[4](111,2) = 1; - this->restriction[4](113,84) = 1; - this->restriction[4](114,39) = 1; - this->restriction[5](5,5) = 1; - this->restriction[5](21,4) = 1; - this->restriction[5](22,21) = 1; - this->restriction[5](23,24) = 1; - this->restriction[5](24,6) = 1; - this->restriction[5](36,1) = 1; - this->restriction[5](37,36) = 1; - this->restriction[5](54,30) = 1; - this->restriction[5](55,57) = 1; - this->restriction[5](57,7) = 1; - this->restriction[5](58,27) = 1; - this->restriction[5](66,0) = 1; - this->restriction[5](67,9) = 1; - this->restriction[5](69,33) = 1; - this->restriction[5](70,66) = 1; - this->restriction[5](72,12) = 1; - this->restriction[5](73,75) = 1; - this->restriction[5](75,2) = 1; - this->restriction[5](76,39) = 1; - this->restriction[5](102,18) = 1; - this->restriction[5](103,48) = 1; - this->restriction[5](105,93) = 1; - this->restriction[5](106,111) = 1; - this->restriction[5](111,3) = 1; - this->restriction[5](112,15) = 1; - this->restriction[5](114,42) = 1; - this->restriction[5](115,84) = 1; - this->restriction[6](6,6) = 1; - this->restriction[6](24,5) = 1; - this->restriction[6](25,24) = 1; - this->restriction[6](27,7) = 1; - this->restriction[6](28,27) = 1; - this->restriction[6](39,2) = 1; - this->restriction[6](40,39) = 1; - this->restriction[6](57,4) = 1; - this->restriction[6](58,21) = 1; - this->restriction[6](60,30) = 1; - this->restriction[6](61,57) = 1; - this->restriction[6](75,1) = 1; - this->restriction[6](76,36) = 1; - this->restriction[6](78,12) = 1; - this->restriction[6](79,75) = 1; - this->restriction[6](84,3) = 1; - this->restriction[6](85,15) = 1; - this->restriction[6](87,42) = 1; - this->restriction[6](88,84) = 1; - this->restriction[6](111,0) = 1; - this->restriction[6](112,9) = 1; - this->restriction[6](114,33) = 1; - this->restriction[6](115,66) = 1; - this->restriction[6](120,18) = 1; - this->restriction[6](121,48) = 1; - this->restriction[6](123,93) = 1; - this->restriction[6](124,111) = 1; - this->restriction[7](7,7) = 1; - this->restriction[7](26,27) = 1; - this->restriction[7](27,6) = 1; - this->restriction[7](30,4) = 1; - this->restriction[7](31,30) = 1; - this->restriction[7](42,3) = 1; - this->restriction[7](43,42) = 1; - this->restriction[7](56,21) = 1; - this->restriction[7](57,5) = 1; - this->restriction[7](59,57) = 1; - this->restriction[7](60,24) = 1; - this->restriction[7](83,15) = 1; - this->restriction[7](84,2) = 1; - this->restriction[7](86,84) = 1; - this->restriction[7](87,39) = 1; - this->restriction[7](93,0) = 1; - this->restriction[7](94,33) = 1; - this->restriction[7](96,18) = 1; - this->restriction[7](97,93) = 1; - this->restriction[7](110,9) = 1; - this->restriction[7](111,1) = 1; - this->restriction[7](113,66) = 1; - this->restriction[7](114,36) = 1; - this->restriction[7](119,48) = 1; - this->restriction[7](120,12) = 1; - this->restriction[7](122,111) = 1; - this->restriction[7](123,75) = 1; - break; - default: - { - // in case we don't - // have the matrices - // (yet), reset them to - // zero size. this does - // not prevent the use - // of this FE, but will - // prevent the use of - // these matrices - for (unsigned int i=0; - i::children_per_cell; - ++i) - this->restriction[i].reinit(0,0); - }; - } - break; - }; - - default: - Assert (false, ExcNotImplemented()); - } - - // finally fill in support points - // on cell and face - initialize_unit_support_points (); - initialize_unit_face_support_points (); -} - - - -template -FiniteElement * -FE_Q::clone() const -{ - return new FE_Q(degree); -} - - - -template -double -FE_Q::shape_value (const unsigned int i, - const Point &p) const -{ - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); - return polynomial_space.compute_value(renumber_inverse[i], p); -} - - -template -double -FE_Q::shape_value_component (const unsigned int i, - const Point &p, - const unsigned int component) const -{ - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); - Assert (component == 0, ExcIndexRange (component, 0, 1)); - return polynomial_space.compute_value(renumber_inverse[i], p); -} - - - -template -Tensor<1,dim> -FE_Q::shape_grad (const unsigned int i, - const Point &p) const -{ - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); - return polynomial_space.compute_grad(renumber_inverse[i], p); -} - - - -template -Tensor<1,dim> -FE_Q::shape_grad_component (const unsigned int i, - const Point &p, - const unsigned int component) const -{ - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); - Assert (component == 0, ExcIndexRange (component, 0, 1)); - return polynomial_space.compute_grad(renumber_inverse[i], p); -} - - - -template -Tensor<2,dim> -FE_Q::shape_grad_grad (const unsigned int i, - const Point &p) const -{ - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); - return polynomial_space.compute_grad_grad(renumber_inverse[i], p); -} - - - -template -Tensor<2,dim> -FE_Q::shape_grad_grad_component (const unsigned int i, - const Point &p, - const unsigned int component) const -{ - Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); - Assert (component == 0, ExcIndexRange (component, 0, 1)); - return polynomial_space.compute_grad_grad(renumber_inverse[i], p); -} - - -//---------------------------------------------------------------------- -// Auxiliary functions -//---------------------------------------------------------------------- - - - -template -void FE_Q::initialize_unit_support_points () -{ - // number of points: (degree+1)^dim - unsigned int n = degree+1; - for (unsigned int i=1; iunit_support_points.resize(n); - - const double step = 1./degree; - Point p; - - unsigned int k=0; - for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz) - for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy) - for (unsigned int ix=0; ix<=degree; ++ix) - { - p(0) = ix * step; - if (dim>1) - p(1) = iy * step; - if (dim>2) - p(2) = iz * step; - - this->unit_support_points[renumber[k++]] = p; - }; -} - - -#if deal_II_dimension == 1 - -template <> -void FE_Q<1>::initialize_unit_face_support_points () -{ - // no faces in 1d, so nothing to do -} - -#endif - - -template -void FE_Q::initialize_unit_face_support_points () -{ - const unsigned int codim = dim-1; - - // number of points: (degree+1)^codim - unsigned int n = degree+1; - for (unsigned int i=1; iunit_face_support_points.resize(n); - - const double step = 1./degree; - Point p; - - unsigned int k=0; - for (unsigned int iz=0; iz <= ((codim>2) ? degree : 0) ; ++iz) - for (unsigned int iy=0; iy <= ((codim>1) ? degree : 0) ; ++iy) - for (unsigned int ix=0; ix<=degree; ++ix) - { - p(0) = ix * step; - if (codim>1) - p(1) = iy * step; - if (codim>2) - p(2) = iz * step; - - this->unit_face_support_points[face_renumber[k++]] = p; - }; -} - - - -template -std::vector -FE_Q::get_dpo_vector(const unsigned int deg) -{ - std::vector dpo(dim+1, static_cast(1)); - for (unsigned int i=1; i -std::vector -FE_Q::lexicographic_to_hierarchic_numbering (const FiniteElementData &fe_data, - const unsigned int degree) -{ - std::vector renumber (fe_data.dofs_per_cell); - - const unsigned int n = degree+1; - - - if (degree == 0) - { - Assert ((fe_data.dofs_per_vertex == 0) && - ((fe_data.dofs_per_line == 0) || (dim == 1)) && - ((fe_data.dofs_per_quad == 0) || (dim == 2)) && - ((fe_data.dofs_per_hex == 0) || (dim == 3)), - ExcInternalError()); - renumber[0] = 0; - }; - - if (degree > 0) + if (degree > 0) { Assert (fe_data.dofs_per_vertex == 1, ExcInternalError()); for (unsigned int i=0; i::vertices_per_cell; ++i) @@ -1015,46 +449,644 @@ FE_Q::lexicographic_to_hierarchic_numbering (const FiniteElementData & } } - if (GeometryInfo::hexes_per_cell > 0) - for (int i=0; i(GeometryInfo::hexes_per_cell); ++i) - { - unsigned int index = fe_data.first_hex_index; - - for (unsigned int jz = 1; jz::hexes_per_cell > 0) + for (int i=0; i(GeometryInfo::hexes_per_cell); ++i) + { + unsigned int index = fe_data.first_hex_index; + + for (unsigned int jz = 1; jz +std::vector +FE_Q::face_lexicographic_to_hierarchic_numbering (const unsigned int degree) +{ + const FiniteElementData fe_data(FE_Q::get_dpo_vector(degree),1); + return FE_Q::lexicographic_to_hierarchic_numbering (fe_data, degree); +} + + +#if (deal_II_dimension == 1) + +template <> +std::vector +FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int) +{ + return std::vector(); +} + +#endif + + + +template +void +FE_Q::initialize_constraints () +{ + // copy constraint matrices if they + // are defined. otherwise leave them + // at invalid size + if ((dim > 1) && (degree < Matrices::n_constraint_matrices+1)) + { + this->interface_constraints. + TableBase<2,double>::reinit (this->interface_constraints_size()); + this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]); + }; +} + + + +template +void +FE_Q::initialize_embedding () +{ + // copy over embedding matrices if + // they are defined + if ((degree < Matrices::n_embedding_matrices+1) && + (Matrices::embedding[degree-1][0] != 0)) + for (unsigned int c=0; c::children_per_cell; ++c) + { + this->prolongation[c].reinit (this->dofs_per_cell, + this->dofs_per_cell); + this->prolongation[c].fill (Matrices::embedding[degree-1][c]); + }; +} + + + +template +void +FE_Q::initialize_restriction () +{ + + // then fill restriction + // matrices. they are hardcoded for + // the first few elements. in + // contrast to the other matrices, + // these are not stored in the + // files fe_q_[123]d.cc, since they + // contain only a rather small + // number of zeros, and storing + // them element-wise is more + // expensive than just setting the + // nonzero elements as done here + for (unsigned int c=0; c::children_per_cell; ++c) + this->restriction[c].reinit (this->dofs_per_cell, this->dofs_per_cell); + switch (dim) + { + case 1: // 1d + { + switch (degree) + { + case 1: + this->restriction[0](0,0) = 1; + this->restriction[1](1,1) = 1; + break; + case 2: + this->restriction[0](0,0) = 1; + this->restriction[0](2,1) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](1,1) = 1; + break; + case 3: + this->restriction[0](0,0) = 1; + this->restriction[0](2,3) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](3,2) = 1; + break; + case 4: + this->restriction[0](0,0) = 1; + this->restriction[0](2,3) = 1; + this->restriction[0](3,1) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](3,0) = 1; + this->restriction[1](4,3) = 1; + break; + + default: + { + // in case we don't + // have the matrices + // (yet), reset them to + // zero size. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + this->restriction[i].reinit(0,0); + }; + } + break; + }; + + case 2: // 2d + { + switch (degree) + { + case 1: + this->restriction[0](0,0) = 1; + this->restriction[1](1,1) = 1; + this->restriction[2](2,2) = 1; + this->restriction[3](3,3) = 1; + break; + case 2: + this->restriction[0](0,0) = 1; + this->restriction[0](4,1) = 1; + this->restriction[0](7,3) = 1; + this->restriction[0](8,2) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](4,0) = 1; + this->restriction[1](5,2) = 1; + this->restriction[1](8,3) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](5,1) = 1; + this->restriction[2](6,3) = 1; + this->restriction[2](8,0) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](6,2) = 1; + this->restriction[3](7,0) = 1; + this->restriction[3](8,1) = 1; + break; + case 3: + this->restriction[0](0,0) = 1; + this->restriction[0](4,5) = 1; + this->restriction[0](10,11) = 1; + this->restriction[0](12,15) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](5,4) = 1; + this->restriction[1](6,7) = 1; + this->restriction[1](13,14) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](7,6) = 1; + this->restriction[2](9,8) = 1; + this->restriction[2](15,12) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](8,9) = 1; + this->restriction[3](11,10) = 1; + this->restriction[3](14,13) = 1; + break; + case 4: + this->restriction[0](0,0) = 1; + this->restriction[0](4,5) = 1; + this->restriction[0](5,1) = 1; + this->restriction[0](13,14) = 1; + this->restriction[0](14,3) = 1; + this->restriction[0](16,20) = 1; + this->restriction[0](17,8) = 1; + this->restriction[0](19,11) = 1; + this->restriction[0](20,2) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](5,0) = 1; + this->restriction[1](6,5) = 1; + this->restriction[1](7,8) = 1; + this->restriction[1](8,2) = 1; + this->restriction[1](17,14) = 1; + this->restriction[1](18,20) = 1; + this->restriction[1](20,3) = 1; + this->restriction[1](21,11) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](8,1) = 1; + this->restriction[2](9,8) = 1; + this->restriction[2](11,3) = 1; + this->restriction[2](12,11) = 1; + this->restriction[2](20,0) = 1; + this->restriction[2](21,5) = 1; + this->restriction[2](23,14) = 1; + this->restriction[2](24,20) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](10,11) = 1; + this->restriction[3](11,2) = 1; + this->restriction[3](14,0) = 1; + this->restriction[3](15,14) = 1; + this->restriction[3](19,5) = 1; + this->restriction[3](20,1) = 1; + this->restriction[3](22,20) = 1; + this->restriction[3](23,8) = 1; + break; + + default: + { + // in case we don't + // have the matrices + // (yet), reset them to + // zero size. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + this->restriction[i].reinit(0,0); + }; + } + break; + }; + + case 3: // 3d + { + switch (degree) + { + case 1: + this->restriction[0](0,0) = 1; + this->restriction[1](1,1) = 1; + this->restriction[2](2,2) = 1; + this->restriction[3](3,3) = 1; + this->restriction[4](4,4) = 1; + this->restriction[5](5,5) = 1; + this->restriction[6](6,6) = 1; + this->restriction[7](7,7) = 1; + break; + case 2: + this->restriction[0](0,0) = 1; + this->restriction[0](8,1) = 1; + this->restriction[0](11,3) = 1; + this->restriction[0](16,4) = 1; + this->restriction[0](20,2) = 1; + this->restriction[0](22,5) = 1; + this->restriction[0](25,7) = 1; + this->restriction[0](26,6) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](8,0) = 1; + this->restriction[1](9,2) = 1; + this->restriction[1](17,5) = 1; + this->restriction[1](20,3) = 1; + this->restriction[1](22,4) = 1; + this->restriction[1](23,6) = 1; + this->restriction[1](26,7) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](9,1) = 1; + this->restriction[2](10,3) = 1; + this->restriction[2](18,6) = 1; + this->restriction[2](20,0) = 1; + this->restriction[2](23,5) = 1; + this->restriction[2](24,7) = 1; + this->restriction[2](26,4) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](10,2) = 1; + this->restriction[3](11,0) = 1; + this->restriction[3](19,7) = 1; + this->restriction[3](20,1) = 1; + this->restriction[3](24,6) = 1; + this->restriction[3](25,4) = 1; + this->restriction[3](26,5) = 1; + this->restriction[4](4,4) = 1; + this->restriction[4](12,5) = 1; + this->restriction[4](15,7) = 1; + this->restriction[4](16,0) = 1; + this->restriction[4](21,6) = 1; + this->restriction[4](22,1) = 1; + this->restriction[4](25,3) = 1; + this->restriction[4](26,2) = 1; + this->restriction[5](5,5) = 1; + this->restriction[5](12,4) = 1; + this->restriction[5](13,6) = 1; + this->restriction[5](17,1) = 1; + this->restriction[5](21,7) = 1; + this->restriction[5](22,0) = 1; + this->restriction[5](23,2) = 1; + this->restriction[5](26,3) = 1; + this->restriction[6](6,6) = 1; + this->restriction[6](13,5) = 1; + this->restriction[6](14,7) = 1; + this->restriction[6](18,2) = 1; + this->restriction[6](21,4) = 1; + this->restriction[6](23,1) = 1; + this->restriction[6](24,3) = 1; + this->restriction[6](26,0) = 1; + this->restriction[7](7,7) = 1; + this->restriction[7](14,6) = 1; + this->restriction[7](15,4) = 1; + this->restriction[7](19,3) = 1; + this->restriction[7](21,5) = 1; + this->restriction[7](24,2) = 1; + this->restriction[7](25,0) = 1; + this->restriction[7](26,1) = 1; + break; + case 3: + this->restriction[0](0,0) = 1; + this->restriction[0](8,9) = 1; + this->restriction[0](14,15) = 1; + this->restriction[0](24,25) = 1; + this->restriction[0](32,35) = 1; + this->restriction[0](40,43) = 1; + this->restriction[0](52,55) = 1; + this->restriction[0](56,63) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](9,8) = 1; + this->restriction[1](10,11) = 1; + this->restriction[1](26,27) = 1; + this->restriction[1](33,34) = 1; + this->restriction[1](41,42) = 1; + this->restriction[1](44,47) = 1; + this->restriction[1](57,62) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](11,10) = 1; + this->restriction[2](13,12) = 1; + this->restriction[2](28,29) = 1; + this->restriction[2](35,32) = 1; + this->restriction[2](46,45) = 1; + this->restriction[2](49,50) = 1; + this->restriction[2](61,58) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](12,13) = 1; + this->restriction[3](15,14) = 1; + this->restriction[3](30,31) = 1; + this->restriction[3](34,33) = 1; + this->restriction[3](48,51) = 1; + this->restriction[3](54,53) = 1; + this->restriction[3](60,59) = 1; + this->restriction[4](4,4) = 1; + this->restriction[4](16,17) = 1; + this->restriction[4](22,23) = 1; + this->restriction[4](25,24) = 1; + this->restriction[4](36,39) = 1; + this->restriction[4](42,41) = 1; + this->restriction[4](53,54) = 1; + this->restriction[4](58,61) = 1; + this->restriction[5](5,5) = 1; + this->restriction[5](17,16) = 1; + this->restriction[5](18,19) = 1; + this->restriction[5](27,26) = 1; + this->restriction[5](37,38) = 1; + this->restriction[5](43,40) = 1; + this->restriction[5](45,46) = 1; + this->restriction[5](59,60) = 1; + this->restriction[6](6,6) = 1; + this->restriction[6](19,18) = 1; + this->restriction[6](21,20) = 1; + this->restriction[6](29,28) = 1; + this->restriction[6](39,36) = 1; + this->restriction[6](47,44) = 1; + this->restriction[6](51,48) = 1; + this->restriction[6](63,56) = 1; + this->restriction[7](7,7) = 1; + this->restriction[7](20,21) = 1; + this->restriction[7](23,22) = 1; + this->restriction[7](31,30) = 1; + this->restriction[7](38,37) = 1; + this->restriction[7](50,49) = 1; + this->restriction[7](55,52) = 1; + this->restriction[7](62,57) = 1; + break; + case 4: + this->restriction[0](0,0) = 1; + this->restriction[0](8,9) = 1; + this->restriction[0](9,1) = 1; + this->restriction[0](17,18) = 1; + this->restriction[0](18,3) = 1; + this->restriction[0](32,33) = 1; + this->restriction[0](33,4) = 1; + this->restriction[0](44,48) = 1; + this->restriction[0](45,12) = 1; + this->restriction[0](47,15) = 1; + this->restriction[0](48,2) = 1; + this->restriction[0](62,66) = 1; + this->restriction[0](63,36) = 1; + this->restriction[0](65,21) = 1; + this->restriction[0](66,5) = 1; + this->restriction[0](89,93) = 1; + this->restriction[0](90,30) = 1; + this->restriction[0](92,42) = 1; + this->restriction[0](93,7) = 1; + this->restriction[0](98,111) = 1; + this->restriction[0](99,75) = 1; + this->restriction[0](101,57) = 1; + this->restriction[0](102,24) = 1; + this->restriction[0](107,84) = 1; + this->restriction[0](108,39) = 1; + this->restriction[0](110,27) = 1; + this->restriction[0](111,6) = 1; + this->restriction[1](1,1) = 1; + this->restriction[1](9,0) = 1; + this->restriction[1](10,9) = 1; + this->restriction[1](11,12) = 1; + this->restriction[1](12,2) = 1; + this->restriction[1](35,36) = 1; + this->restriction[1](36,5) = 1; + this->restriction[1](45,18) = 1; + this->restriction[1](46,48) = 1; + this->restriction[1](48,3) = 1; + this->restriction[1](49,15) = 1; + this->restriction[1](63,33) = 1; + this->restriction[1](64,66) = 1; + this->restriction[1](66,4) = 1; + this->restriction[1](67,21) = 1; + this->restriction[1](71,75) = 1; + this->restriction[1](72,24) = 1; + this->restriction[1](74,39) = 1; + this->restriction[1](75,6) = 1; + this->restriction[1](99,93) = 1; + this->restriction[1](100,111) = 1; + this->restriction[1](102,30) = 1; + this->restriction[1](103,57) = 1; + this->restriction[1](108,42) = 1; + this->restriction[1](109,84) = 1; + this->restriction[1](111,7) = 1; + this->restriction[1](112,27) = 1; + this->restriction[2](2,2) = 1; + this->restriction[2](12,1) = 1; + this->restriction[2](13,12) = 1; + this->restriction[2](15,3) = 1; + this->restriction[2](16,15) = 1; + this->restriction[2](38,39) = 1; + this->restriction[2](39,6) = 1; + this->restriction[2](48,0) = 1; + this->restriction[2](49,9) = 1; + this->restriction[2](51,18) = 1; + this->restriction[2](52,48) = 1; + this->restriction[2](74,36) = 1; + this->restriction[2](75,5) = 1; + this->restriction[2](77,75) = 1; + this->restriction[2](78,24) = 1; + this->restriction[2](81,42) = 1; + this->restriction[2](82,84) = 1; + this->restriction[2](84,7) = 1; + this->restriction[2](85,27) = 1; + this->restriction[2](108,33) = 1; + this->restriction[2](109,66) = 1; + this->restriction[2](111,4) = 1; + this->restriction[2](112,21) = 1; + this->restriction[2](117,93) = 1; + this->restriction[2](118,111) = 1; + this->restriction[2](120,30) = 1; + this->restriction[2](121,57) = 1; + this->restriction[3](3,3) = 1; + this->restriction[3](14,15) = 1; + this->restriction[3](15,2) = 1; + this->restriction[3](18,0) = 1; + this->restriction[3](19,18) = 1; + this->restriction[3](41,42) = 1; + this->restriction[3](42,7) = 1; + this->restriction[3](47,9) = 1; + this->restriction[3](48,1) = 1; + this->restriction[3](50,48) = 1; + this->restriction[3](51,12) = 1; + this->restriction[3](80,84) = 1; + this->restriction[3](81,39) = 1; + this->restriction[3](83,27) = 1; + this->restriction[3](84,6) = 1; + this->restriction[3](92,33) = 1; + this->restriction[3](93,4) = 1; + this->restriction[3](95,93) = 1; + this->restriction[3](96,30) = 1; + this->restriction[3](107,66) = 1; + this->restriction[3](108,36) = 1; + this->restriction[3](110,21) = 1; + this->restriction[3](111,5) = 1; + this->restriction[3](116,111) = 1; + this->restriction[3](117,75) = 1; + this->restriction[3](119,57) = 1; + this->restriction[3](120,24) = 1; + this->restriction[4](4,4) = 1; + this->restriction[4](20,21) = 1; + this->restriction[4](21,5) = 1; + this->restriction[4](29,30) = 1; + this->restriction[4](30,7) = 1; + this->restriction[4](33,0) = 1; + this->restriction[4](34,33) = 1; + this->restriction[4](53,57) = 1; + this->restriction[4](54,24) = 1; + this->restriction[4](56,27) = 1; + this->restriction[4](57,6) = 1; + this->restriction[4](65,9) = 1; + this->restriction[4](66,1) = 1; + this->restriction[4](68,66) = 1; + this->restriction[4](69,36) = 1; + this->restriction[4](90,18) = 1; + this->restriction[4](91,93) = 1; + this->restriction[4](93,3) = 1; + this->restriction[4](94,42) = 1; + this->restriction[4](101,48) = 1; + this->restriction[4](102,12) = 1; + this->restriction[4](104,111) = 1; + this->restriction[4](105,75) = 1; + this->restriction[4](110,15) = 1; + this->restriction[4](111,2) = 1; + this->restriction[4](113,84) = 1; + this->restriction[4](114,39) = 1; + this->restriction[5](5,5) = 1; + this->restriction[5](21,4) = 1; + this->restriction[5](22,21) = 1; + this->restriction[5](23,24) = 1; + this->restriction[5](24,6) = 1; + this->restriction[5](36,1) = 1; + this->restriction[5](37,36) = 1; + this->restriction[5](54,30) = 1; + this->restriction[5](55,57) = 1; + this->restriction[5](57,7) = 1; + this->restriction[5](58,27) = 1; + this->restriction[5](66,0) = 1; + this->restriction[5](67,9) = 1; + this->restriction[5](69,33) = 1; + this->restriction[5](70,66) = 1; + this->restriction[5](72,12) = 1; + this->restriction[5](73,75) = 1; + this->restriction[5](75,2) = 1; + this->restriction[5](76,39) = 1; + this->restriction[5](102,18) = 1; + this->restriction[5](103,48) = 1; + this->restriction[5](105,93) = 1; + this->restriction[5](106,111) = 1; + this->restriction[5](111,3) = 1; + this->restriction[5](112,15) = 1; + this->restriction[5](114,42) = 1; + this->restriction[5](115,84) = 1; + this->restriction[6](6,6) = 1; + this->restriction[6](24,5) = 1; + this->restriction[6](25,24) = 1; + this->restriction[6](27,7) = 1; + this->restriction[6](28,27) = 1; + this->restriction[6](39,2) = 1; + this->restriction[6](40,39) = 1; + this->restriction[6](57,4) = 1; + this->restriction[6](58,21) = 1; + this->restriction[6](60,30) = 1; + this->restriction[6](61,57) = 1; + this->restriction[6](75,1) = 1; + this->restriction[6](76,36) = 1; + this->restriction[6](78,12) = 1; + this->restriction[6](79,75) = 1; + this->restriction[6](84,3) = 1; + this->restriction[6](85,15) = 1; + this->restriction[6](87,42) = 1; + this->restriction[6](88,84) = 1; + this->restriction[6](111,0) = 1; + this->restriction[6](112,9) = 1; + this->restriction[6](114,33) = 1; + this->restriction[6](115,66) = 1; + this->restriction[6](120,18) = 1; + this->restriction[6](121,48) = 1; + this->restriction[6](123,93) = 1; + this->restriction[6](124,111) = 1; + this->restriction[7](7,7) = 1; + this->restriction[7](26,27) = 1; + this->restriction[7](27,6) = 1; + this->restriction[7](30,4) = 1; + this->restriction[7](31,30) = 1; + this->restriction[7](42,3) = 1; + this->restriction[7](43,42) = 1; + this->restriction[7](56,21) = 1; + this->restriction[7](57,5) = 1; + this->restriction[7](59,57) = 1; + this->restriction[7](60,24) = 1; + this->restriction[7](83,15) = 1; + this->restriction[7](84,2) = 1; + this->restriction[7](86,84) = 1; + this->restriction[7](87,39) = 1; + this->restriction[7](93,0) = 1; + this->restriction[7](94,33) = 1; + this->restriction[7](96,18) = 1; + this->restriction[7](97,93) = 1; + this->restriction[7](110,9) = 1; + this->restriction[7](111,1) = 1; + this->restriction[7](113,66) = 1; + this->restriction[7](114,36) = 1; + this->restriction[7](119,48) = 1; + this->restriction[7](120,12) = 1; + this->restriction[7](122,111) = 1; + this->restriction[7](123,75) = 1; + break; + default: + { + // in case we don't + // have the matrices + // (yet), reset them to + // zero size. this does + // not prevent the use + // of this FE, but will + // prevent the use of + // these matrices + for (unsigned int i=0; + i::children_per_cell; + ++i) + this->restriction[i].reinit(0,0); + }; + } + break; + }; + + default: + Assert (false, ExcNotImplemented()); } - - return renumber; -} - - - -template -std::vector -FE_Q::face_lexicographic_to_hierarchic_numbering (const unsigned int degree) -{ - const FiniteElementData fe_data(FE_Q::get_dpo_vector(degree),1); - return FE_Q::lexicographic_to_hierarchic_numbering (fe_data, degree); -} - - -#if (deal_II_dimension == 1) - -template <> -std::vector -FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int) -{ - return std::vector(); } -#endif template -- 2.39.5