From: wolf Date: Mon, 30 Sep 2002 14:07:58 +0000 (+0000) Subject: No longer preset the correct sizes of matrices in the FE base class, but force those... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0f6f10d7410f0c724e310f9f0b1f5605fabff390;p=dealii-svn.git No longer preset the correct sizes of matrices in the FE base class, but force those classes that implement these matrices to set their size. This way, classes that do _not_ implement them, cannot forget to reset their size to zero, as has happened in one of the discontinuous elements. git-svn-id: https://svn.dealii.org/trunk@6557 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_base.h b/deal.II/deal.II/include/fe/fe_base.h index 271eb47321..297e0846ac 100644 --- a/deal.II/deal.II/include/fe/fe_base.h +++ b/deal.II/deal.II/include/fe/fe_base.h @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -1106,6 +1107,21 @@ class FiniteElementBase : public Subscriptor, */ FullMatrix interface_constraints; + /** + * Return the size of interface + * constraint matrices. Since + * this is needed in every + * derived finite element class + * when initializing their size, + * it is placed into this + * function, to avoid having to + * recomputed the + * dimension-dependent size of + * these matrices each time. + */ + TableIndices<2> + interface_constraints_size () const; + /** * Store what * @p{system_to_component_index} diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index cc08011d62..9c59f139e5 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -27,6 +27,9 @@ #include +//TODO[GK]: in all FEs, somehow make the handling of first_cell more intuitive. presently, this flag is cleared at the end of the fill_fe_*_values function, but this leads to confusion, if for example in the FESystems class, we ask for the present values of the update flags, using update_current, because the element thinks it is already on the second cell and the FESystem does not want to do some actions it is supposed to do. thus, in the respective FESystem function, we store these update_values once before calling the fill_* functions of the subelements, and reuse the stored values. if we would request their values only after initializing the base element (i.e. requesting them at the point _where we actually need them_!), this would fail, and we would not copy the shape_values, for example. +//TODO[GK]: the solution would probably to not set and clear the first_cell flag in the FEs themselves, but to have a function which sets or clears this flag, and call it from the FEValues object, which should always know what the present state is. this way, we could be more sure that this flag is actually consistent with the present state, and among the various elements that are involved in a FESystem. + // if necessary try to work around a bug in the IBM xlC compiler #ifdef XLC_WORK_AROUND_STD_BUG using namespace std; @@ -151,43 +154,9 @@ FiniteElementBase (const FiniteElementData &fe_data, ExcInternalError()); }; - for (unsigned int i=0; i::children_per_cell; ++i) - { - restriction[i].reinit (this->dofs_per_cell, this->dofs_per_cell); - prolongation[i].reinit (this->dofs_per_cell, this->dofs_per_cell); - }; - - // first set sizes of some - // matrices. they will be filled by - // derived classes, or re-set to - // zero size of not defined - switch (dim) - { - case 1: - Assert ((interface_constraints.m() == 0) && - (interface_constraints.n() == 0), - ExcInternalError()); - break; - - case 2: - interface_constraints.reinit (this->dofs_per_vertex - +2*this->dofs_per_line, - this->dofs_per_face); - break; - - case 3: - interface_constraints.reinit (5*this->dofs_per_vertex + - 12*this->dofs_per_line + - 4*this->dofs_per_quad, - this->dofs_per_face); - break; - - default: - Assert (false, ExcNotImplemented()); - }; - - // this is the default way, if - // there is only one component; if + // initialize some tables in the + // default way, i.e. if there is + // only one (vector-)component; if // there are several, then the // constructor of the derived class // needs to overwrite these arrays @@ -313,6 +282,32 @@ FiniteElementBase::constraints () const +template +TableIndices<2> +FiniteElementBase::interface_constraints_size () const +{ + switch (dim) + { + case 1: + return TableIndices<2> (0U, 0U); + case 2: + return TableIndices<2> (this->dofs_per_vertex + + 2*this->dofs_per_line, + this->dofs_per_face); + case 3: + return TableIndices<2> (5*this->dofs_per_vertex + + 12*this->dofs_per_line + + 4*this->dofs_per_quad, + this->dofs_per_face); + default: + Assert (false, ExcNotImplemented()); + }; + return TableIndices<2> (-1, -1); +}; + + + + template bool FiniteElementBase::operator == (const FiniteElementBase &f) const { diff --git a/deal.II/deal.II/source/fe/fe_dgp.cc b/deal.II/deal.II/source/fe/fe_dgp.cc index b8331e0993..7018ccf84c 100644 --- a/deal.II/deal.II/source/fe/fe_dgp.cc +++ b/deal.II/deal.II/source/fe/fe_dgp.cc @@ -37,7 +37,11 @@ FE_DGP::FE_DGP (const unsigned int degree) if ((degree < Matrices::n_embedding_matrices) && (Matrices::embedding[degree][0] != 0)) for (unsigned int c=0; c::children_per_cell; ++c) - this->prolongation[c].fill (Matrices::embedding[degree][c]); + { + this->prolongation[c].reinit (this->dofs_per_cell, + this->dofs_per_cell); + this->prolongation[c].fill (Matrices::embedding[degree][c]); + } else for (unsigned int i=0; i::children_per_cell;++i) this->prolongation[i].reinit(0,0); @@ -50,6 +54,7 @@ FE_DGP::FE_DGP (const unsigned int degree) // // if it were, then the following // snippet would be the right code +// this->restriction[i].reinit (this->dofs_per_cell, this->dofs_per_cell); // if ((degree < Matrices::n_projection_matrices) && // (Matrices::projection_matrices[degree] != 0)) // { @@ -59,10 +64,6 @@ FE_DGP::FE_DGP (const unsigned int degree) // // matrix undefined, set size to zero // for (unsigned int i=0;i::children_per_cell;++i) // restriction[i].reinit(0, 0); - // since not implemented, set to - // "empty" - for (unsigned int i=0;i::children_per_cell;++i) - restriction[i].reinit(0, 0); // note further, that these // elements have neither support diff --git a/deal.II/deal.II/source/fe/fe_dgq.cc b/deal.II/deal.II/source/fe/fe_dgq.cc index 6ddd18b930..9f7f38eb8a 100644 --- a/deal.II/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/deal.II/source/fe/fe_dgq.cc @@ -57,14 +57,21 @@ FE_DGQ::FE_DGQ (const unsigned int degree) if ((degree < Matrices::n_embedding_matrices) && (Matrices::embedding[degree] != 0)) { + for (unsigned int i=0; i::children_per_cell; ++i) + this->prolongation[i].reinit (this->dofs_per_cell, + this->dofs_per_cell); this->prolongation[0].fill (Matrices::embedding[degree]); switch (dim) { case 1: + { this->prolongation[1].fill_permutation (this->prolongation[0], right, right); break; + }; + case 2: + { this->prolongation[1].fill_permutation (this->prolongation[0], right, right); this->prolongation[2].fill_permutation (this->prolongation[1], @@ -72,7 +79,10 @@ FE_DGQ::FE_DGQ (const unsigned int degree) this->prolongation[3].fill_permutation (this->prolongation[2], right, right); break; + }; + case 3: + { this->prolongation[1].fill_permutation (this->prolongation[0], right, right); this->prolongation[5].fill_permutation (this->prolongation[1], @@ -88,14 +98,16 @@ FE_DGQ::FE_DGQ (const unsigned int degree) this->prolongation[2].fill_permutation (this->prolongation[6], top, top); break; + }; + default: Assert (false, ExcNotImplemented()); } } else - // matrix undefined, set size to zero - for (unsigned int i=0;i::children_per_cell;++i) - this->prolongation[i].reinit(0, 0); + // matrix undefined, leave matrix + // at size zero + ; // same as above: copy over matrix // from predefined values and @@ -103,14 +115,21 @@ FE_DGQ::FE_DGQ (const unsigned int degree) if ((degree < Matrices::n_projection_matrices) && (Matrices::projection_matrices[degree] != 0)) { + for (unsigned int i=0; i::children_per_cell; ++i) + this->restriction[i].reinit (this->dofs_per_cell, + this->dofs_per_cell); this->restriction[0].fill (Matrices::projection_matrices[degree]); switch (dim) { case 1: + { this->restriction[1].fill_permutation (this->restriction[0], right, right); break; + }; + case 2: + { this->restriction[1].fill_permutation (this->restriction[0], right, right); this->restriction[2].fill_permutation (this->restriction[1], @@ -118,7 +137,10 @@ FE_DGQ::FE_DGQ (const unsigned int degree) this->restriction[3].fill_permutation (this->restriction[2], right, right); break; + }; + case 3: + { this->restriction[1].fill_permutation (this->restriction[0], right, right); this->restriction[5].fill_permutation (this->restriction[1], @@ -134,15 +156,16 @@ FE_DGQ::FE_DGQ (const unsigned int degree) this->restriction[2].fill_permutation (this->restriction[6], top, top); break; + }; + default: Assert (false, ExcNotImplemented()); } } else - // matrix undefined, set size to zero - for (unsigned int i=0;i::children_per_cell;++i) - this->restriction[i].reinit(0, 0); - + // matrix undefined, leave matrix + // at size zero + ; // finally fill in support points if (degree == 0) diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index 22ada26dc7..079b01e0f1 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -12,8 +12,7 @@ //---------------------------------------------------------------- #include -#include -#include +#include #include #include #include @@ -36,12 +35,14 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) Assert (dim >= 2, ExcNotUsefulInThisDimension()); // copy constraint matrices if they - // are defined. otherwise set them - // to invalid size + // are defined. otherwise leave them + // at zero size if (degreeinterface_constraints.fill (Matrices::constraint_matrices[degree-1]); - else - this->interface_constraints.reinit(0,0); + { + this->interface_constraints. + TableBase<2,double>::reinit (this->interface_constraints_size()); + this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]); + }; // next copy over embedding // matrices if they are defined @@ -50,6 +51,8 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) for (unsigned int c=0; c::children_per_cell; ++c) { // copy + this->prolongation[c].reinit (this->dofs_per_cell, + this->dofs_per_cell); this->prolongation[c].fill (Matrices::embedding[degree-1][c]); // and make sure that the row // sum is 0.5 (for usual @@ -68,10 +71,7 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) Assert (std::fabs(sum-.5) < 1e-14, ExcInternalError()); }; - } - else - for (unsigned int i=0; i::children_per_cell;++i) - this->prolongation[i].reinit(0,0); + }; // then fill restriction // matrices. they are hardcoded for @@ -165,6 +165,10 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) // always, in the // canonical direction // of lines + for (unsigned int c=0; c::children_per_cell; ++c) + this->restriction[c].reinit (this->dofs_per_cell, + this->dofs_per_cell); + this->restriction[0](0,0) = 2.; this->restriction[1](1,1) = 2.; this->restriction[3](2,2) = 2.; @@ -177,17 +181,13 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) { // in case we don't // have the matrices - // (yet), set them to - // impossible - // values. this does - // not prevent the use - // of this FE, but will + // (yet), leave them + // empty. 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; }; }; @@ -206,6 +206,9 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) // cell to get at the // values of each of // the 12 lines + for (unsigned int c=0; c::children_per_cell; ++c) + this->restriction[c].reinit (this->dofs_per_cell, + this->dofs_per_cell); this->restriction[0](0,0) = 2.; this->restriction[0](3,3) = 2.; this->restriction[1](1,1) = 2.; @@ -228,17 +231,13 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) { // in case we don't // have the matrices - // (yet), set them to - // impossible - // values. this does - // not prevent the use - // of this FE, but will + // (yet), leave them + // empty. 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; }; }; @@ -253,6 +252,14 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) // on cell and face initialize_unit_support_points (); initialize_unit_face_support_points (); + + // then make + // system_to_component_table + // invalid, since this has no + // meaning for the present element + std::vector > tmp1, tmp2; + this->system_to_component_table.swap (tmp1); + this->face_system_to_component_table.swap (tmp2); }; diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 845629644e..450063a9e5 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -49,22 +49,25 @@ FE_Q::FE_Q (const unsigned int degree) renumber_inverse[renumber[i]]=i; // copy constraint matrices if they - // are defined. otherwise set them - // to invalid size + // are defined. otherwise leave them + // at invalid size if ((dim>1) && (degreeinterface_constraints.fill (Matrices::constraint_matrices[degree-1]); - else - this->interface_constraints.reinit(0,0); + { + this->interface_constraints. + TableBase<2,double>::reinit (this->interface_constraints_size()); + this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]); + }; // 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].fill (Matrices::embedding[degree-1][c]); - else - for (unsigned int i=0; i::children_per_cell;++i) - this->prolongation[i].reinit(0,0); + { + 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 @@ -77,548 +80,547 @@ FE_Q::FE_Q (const unsigned int degree) // 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; + { + 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), set them to - // impossible - // values. 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: + { + // 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), set them to - // impossible - // values. 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; - }; + { + 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), set them to - // impossible - // values. 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; - }; + { + 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()); + Assert (false, ExcNotImplemented()); } // finally fill in support points diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index a890debab6..0579ba33a5 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -1215,6 +1215,9 @@ FESystem::build_face_tables() template void FESystem::build_interface_constraints () { + this->interface_constraints. + TableBase<2,double>::reinit (this->interface_constraints_size()); + // the layout of the constraints // matrix is described in the // FiniteElement class. you may @@ -1482,15 +1485,17 @@ void FESystem::initialize () do_prolongation = false; } - // if we encountered void matrices, - // disable them for the composite - // element as well - if (!do_restriction) + // if we did not encounter void + // matrices, initialize the + // respective matrix sizes + if (do_restriction) for (unsigned int i=0;i::children_per_cell;++i) - this->restriction[i].reinit(0,0); - if (!do_prolongation) + this->restriction[i].reinit(this->dofs_per_cell, + this->dofs_per_cell); + if (do_prolongation) for (unsigned int i=0;i::children_per_cell;++i) - this->prolongation[i].reinit(0,0); + this->prolongation[i].reinit(this->dofs_per_cell, + this->dofs_per_cell); // distribute the matrices of the // base finite elements to the