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<Matrices::n_constraint_matrices+1))
+
+ // copy constraint and embedding
+ // matrices if they are
+ // defined. otherwise leave them at
+ // invalid size
+ initialize_constraints ();
+ initialize_embedding ();
+ initialize_restriction ();
+
+ // finally fill in support points
+ // on cell and face
+ initialize_unit_support_points ();
+ initialize_unit_face_support_points ();
+}
+
+
+
+template <int dim>
+FiniteElement<dim> *
+FE_Q<dim>::clone() const
+{
+ return new FE_Q<dim>(degree);
+}
+
+
+
+template <int dim>
+double
+FE_Q<dim>::shape_value (const unsigned int i,
+ const Point<dim> &p) const
+{
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+ return polynomial_space.compute_value(renumber_inverse[i], p);
+}
+
+
+template <int dim>
+double
+FE_Q<dim>::shape_value_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<this->dofs_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 <int dim>
+Tensor<1,dim>
+FE_Q<dim>::shape_grad (const unsigned int i,
+ const Point<dim> &p) const
+{
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+ return polynomial_space.compute_grad(renumber_inverse[i], p);
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+FE_Q<dim>::shape_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<this->dofs_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 <int dim>
+Tensor<2,dim>
+FE_Q<dim>::shape_grad_grad (const unsigned int i,
+ const Point<dim> &p) const
+{
+ Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+ return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+FE_Q<dim>::shape_grad_grad_component (const unsigned int i,
+ const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (i<this->dofs_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 <int dim>
+void FE_Q<dim>::initialize_unit_support_points ()
+{
+ // number of points: (degree+1)^dim
+ unsigned int n = degree+1;
+ for (unsigned int i=1; i<dim; ++i)
+ n *= degree+1;
+
+ this->unit_support_points.resize(n);
+
+ const double step = 1./degree;
+ Point<dim> 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 <int dim>
+void FE_Q<dim>::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; i<codim; ++i)
+ n *= degree+1;
+
+ this->unit_face_support_points.resize(n);
+
+ const double step = 1./degree;
+ Point<codim> 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 <int dim>
+std::vector<unsigned int>
+FE_Q<dim>::get_dpo_vector(const unsigned int deg)
+{
+ std::vector<unsigned int> dpo(dim+1, static_cast<unsigned int>(1));
+ for (unsigned int i=1; i<dpo.size(); ++i)
+ dpo[i]=dpo[i-1]*(deg-1);
+ return dpo;
+}
+
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_Q<dim>::lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
+ const unsigned int degree)
+{
+ std::vector<unsigned int> 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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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 <int dim>
-FiniteElement<dim> *
-FE_Q<dim>::clone() const
-{
- return new FE_Q<dim>(degree);
-}
-
-
-
-template <int dim>
-double
-FE_Q<dim>::shape_value (const unsigned int i,
- const Point<dim> &p) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
- return polynomial_space.compute_value(renumber_inverse[i], p);
-}
-
-
-template <int dim>
-double
-FE_Q<dim>::shape_value_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_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 <int dim>
-Tensor<1,dim>
-FE_Q<dim>::shape_grad (const unsigned int i,
- const Point<dim> &p) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
- return polynomial_space.compute_grad(renumber_inverse[i], p);
-}
-
-
-
-template <int dim>
-Tensor<1,dim>
-FE_Q<dim>::shape_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_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 <int dim>
-Tensor<2,dim>
-FE_Q<dim>::shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const
-{
- Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
- return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
-}
-
-
-
-template <int dim>
-Tensor<2,dim>
-FE_Q<dim>::shape_grad_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
-{
- Assert (i<this->dofs_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 <int dim>
-void FE_Q<dim>::initialize_unit_support_points ()
-{
- // number of points: (degree+1)^dim
- unsigned int n = degree+1;
- for (unsigned int i=1; i<dim; ++i)
- n *= degree+1;
-
- this->unit_support_points.resize(n);
-
- const double step = 1./degree;
- Point<dim> 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 <int dim>
-void FE_Q<dim>::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; i<codim; ++i)
- n *= degree+1;
-
- this->unit_face_support_points.resize(n);
-
- const double step = 1./degree;
- Point<codim> 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 <int dim>
-std::vector<unsigned int>
-FE_Q<dim>::get_dpo_vector(const unsigned int deg)
-{
- std::vector<unsigned int> dpo(dim+1, static_cast<unsigned int>(1));
- for (unsigned int i=1; i<dpo.size(); ++i)
- dpo[i]=dpo[i-1]*(deg-1);
- return dpo;
-}
-
-
-
-template <int dim>
-std::vector<unsigned int>
-FE_Q<dim>::lexicographic_to_hierarchic_numbering (const FiniteElementData<dim> &fe_data,
- const unsigned int degree)
-{
- std::vector<unsigned int> 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<GeometryInfo<dim>::vertices_per_cell; ++i)
}
}
- if (GeometryInfo<dim>::hexes_per_cell > 0)
- for (int i=0; i<static_cast<signed int>(GeometryInfo<dim>::hexes_per_cell); ++i)
- {
- unsigned int index = fe_data.first_hex_index;
-
- for (unsigned int jz = 1; jz<degree; jz++)
- for (unsigned int jy = 1; jy<degree; jy++)
- for (unsigned int jx = 1; jx<degree; jx++)
- {
- const unsigned int tensorindex = jx + jy*n + jz*n*n;
- Assert (tensorindex<renumber.size(), ExcInternalError());
- renumber[tensorindex]=index++;
- }
- }
+ if (GeometryInfo<dim>::hexes_per_cell > 0)
+ for (int i=0; i<static_cast<signed int>(GeometryInfo<dim>::hexes_per_cell); ++i)
+ {
+ unsigned int index = fe_data.first_hex_index;
+
+ for (unsigned int jz = 1; jz<degree; jz++)
+ for (unsigned int jy = 1; jy<degree; jy++)
+ for (unsigned int jx = 1; jx<degree; jx++)
+ {
+ const unsigned int tensorindex = jx + jy*n + jz*n*n;
+ Assert (tensorindex<renumber.size(), ExcInternalError());
+ renumber[tensorindex]=index++;
+ }
+ }
+ }
+
+ return renumber;
+}
+
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_Q<dim>::face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
+{
+ const FiniteElementData<dim-1> fe_data(FE_Q<dim-1>::get_dpo_vector(degree),1);
+ return FE_Q<dim-1>::lexicographic_to_hierarchic_numbering (fe_data, degree);
+}
+
+
+#if (deal_II_dimension == 1)
+
+template <>
+std::vector<unsigned int>
+FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int)
+{
+ return std::vector<unsigned int>();
+}
+
+#endif
+
+
+
+template <int dim>
+void
+FE_Q<dim>::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 <int dim>
+void
+FE_Q<dim>::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<GeometryInfo<dim>::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 <int dim>
+void
+FE_Q<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ this->restriction[i].reinit(0,0);
+ };
+ }
+ break;
+ };
+
+ default:
+ Assert (false, ExcNotImplemented());
}
-
- return renumber;
-}
-
-
-
-template <int dim>
-std::vector<unsigned int>
-FE_Q<dim>::face_lexicographic_to_hierarchic_numbering (const unsigned int degree)
-{
- const FiniteElementData<dim-1> fe_data(FE_Q<dim-1>::get_dpo_vector(degree),1);
- return FE_Q<dim-1>::lexicographic_to_hierarchic_numbering (fe_data, degree);
-}
-
-
-#if (deal_II_dimension == 1)
-
-template <>
-std::vector<unsigned int>
-FE_Q<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int)
-{
- return std::vector<unsigned int>();
}
-#endif
template <int dim>
projects / dealii-svn.git / commitdiff