From 92585bf5351555f1146d4904df2f00a01e02eb04 Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 23 Apr 2003 23:56:24 +0000 Subject: [PATCH] Unify the way the various finite elements are initialized in their constructors, by having functions of the same name doing the various stuff that is usually done in the constructor. This should make it significantly simpler in the future to start with a new element, since the various steps that need to be taken are now a little clearer. git-svn-id: https://svn.dealii.org/trunk@7464 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_nedelec.h | 22 + .../deal.II/include/fe/fe_q_hierarchical.h | 67 +- deal.II/deal.II/source/fe/fe_nedelec.cc | 450 ++++++------ deal.II/deal.II/source/fe/fe_q.cc | 13 +- .../deal.II/source/fe/fe_q_hierarchical.cc | 656 ++++++++++-------- 5 files changed, 670 insertions(+), 538 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_nedelec.h b/deal.II/deal.II/include/fe/fe_nedelec.h index a7553e53be..e1d1472ca8 100644 --- a/deal.II/deal.II/include/fe/fe_nedelec.h +++ b/deal.II/deal.II/include/fe/fe_nedelec.h @@ -447,6 +447,28 @@ class FE_Nedelec : public FiniteElement */ static std::vector get_dpo_vector(const unsigned int degree); + + /** + * Initialize the hanging node + * constraints matrices. Called + * from the constructor. + */ + void initialize_constraints (); + + /** + * Initialize the embedding + * matrices. Called from the + * constructor. + */ + void initialize_embedding (); + + /** + * Initialize the restriction + * matrices. Called from the + * constructor. + */ + void initialize_restriction (); + /** * Initialize the * @p{unit_support_points} field diff --git a/deal.II/deal.II/include/fe/fe_q_hierarchical.h b/deal.II/deal.II/include/fe/fe_q_hierarchical.h index 54d59b934b..816bb73f68 100644 --- a/deal.II/deal.II/include/fe/fe_q_hierarchical.h +++ b/deal.II/deal.II/include/fe/fe_q_hierarchical.h @@ -560,10 +560,9 @@ class FE_Q_Hierarchical : public FiniteElement * proper design. */ static - void + std::vector lexicographic_to_hierarchic_numbering (const FiniteElementData &fe_data, - const unsigned int degree, - std::vector &numbering); + const unsigned int degree); /** * This is an analogon to the @@ -571,10 +570,32 @@ class FE_Q_Hierarchical : public FiniteElement * on faces. */ static - void - face_lexicographic_to_hierarchic_numbering (const unsigned int degree, - std::vector &numbering); + std::vector + face_lexicographic_to_hierarchic_numbering (const unsigned int degree); + /** + * Initialize two auxiliary + * fields that will be used in + * setting up the various + * matrices in the constructor. + */ + void build_dofs_cell (std::vector > &dofs_cell, + std::vector > &dofs_subcell) const; + + /** + * Initialize the hanging node + * constraints matrices. Called + * from the constructor. + */ + void initialize_constraints (const std::vector > &dofs_subcell); + + /** + * Initialize the embedding + * matrices. Called from the + * constructor. + */ + void initialize_embedding_and_restriction (const std::vector > &dofs_cell, + const std::vector > &dofs_subcell); /** * Initialize the @@ -672,7 +693,7 @@ class FE_Q_Hierarchical : public FiniteElement * Mapping from lexicographic to * shape function numbering. */ - std::vector renumber; + const std::vector renumber; /** * Inverse renumber @@ -680,32 +701,13 @@ class FE_Q_Hierarchical : public FiniteElement * shape function numbering to * lexicographic numbering. */ - std::vector renumber_inverse; + const std::vector renumber_inverse; /** * Mapping from lexicographic to * shape function numbering on first face. */ - std::vector face_renumber; - - /** - * The matrix @p{dofs_cell} contains the - * values of the linear functionals of - * the master 1d cell applied to the - * shape functions of the two 1d subcells. - * The matrix @p{dofs_subcell} constains - * the values of the linear functionals - * on each 1d subcell applied to the - * shape functions on the master 1d - * subcell. - * We use @p{dofs_cell} and - * @p{dofs_subcell} to compute the - * @p{prolongation}, @p{restriction} and - * @p{interface_constraints} matrices - * for all dimensions. - */ - std::vector > dofs_cell; - std::vector > dofs_subcell; + const std::vector face_renumber; /** * Pointer to the tensor @@ -782,8 +784,11 @@ class FE_Q_Hierarchical : public FiniteElement /* -------------- declaration of explicit specializations ------------- */ -template <> void FE_Q_Hierarchical<1>::initialize_unit_face_support_points (); -template <> void FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int, - std::vector&); +template <> +void FE_Q_Hierarchical<1>::initialize_unit_face_support_points (); + +template <> +std::vector +FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int); #endif diff --git a/deal.II/deal.II/source/fe/fe_nedelec.cc b/deal.II/deal.II/source/fe/fe_nedelec.cc index f93bc4e22f..37bf020088 100644 --- a/deal.II/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/deal.II/source/fe/fe_nedelec.cc @@ -34,219 +34,13 @@ FE_Nedelec::FE_Nedelec (const unsigned int degree) { Assert (dim >= 2, ExcNotUsefulInThisDimension()); - // copy constraint matrices if they - // are defined. otherwise leave them - // at zero size - if (degreeinterface_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) - { - // 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 - // elements, the row sum must - // be 1, but here the shape - // function is multiplied by - // the inverse of the - // Jacobian, which introduces - // a factor of 1/2 when going - // from mother to child) - for (unsigned int row=0; rowdofs_per_cell; ++row) - { - double sum = 0; - for (unsigned int col=0; coldofs_per_cell; ++col) - sum += this->prolongation[c](row,col); - Assert (std::fabs(sum-.5) < 1e-14, - ExcInternalError()); - }; - }; - - // then fill restriction - // matrices. they are hardcoded for - // the first few elements - switch (dim) - { - case 2: // 2d - { - switch (degree) - { - case 1: - { - // this is a strange - // element, since it is - // both additive and - // then it is also - // not. ideally, we - // would like to have - // the value of the - // shape function on - // the coarse line to - // be the mean value of - // that on the two - // child ones. thus, - // one should make it - // additive. however, - // additivity only - // works if an element - // does not have any - // continuity - // requirements, since - // otherwise degrees of - // freedom are shared - // between adjacent - // elements, and when - // we make the element - // additive, that would - // mean that we end up - // adding up - // contributions not - // only from the child - // cells of this cell, - // but also from the - // child cells of the - // neighbor, and since - // we cannot know - // whether there even - // exists a neighbor we - // cannot simply make - // the element - // additive. - // - // so, until someone - // comes along with a - // better alternative, - // we do the following: - // make the element - // non-additive, and - // simply pick the - // value of one of the - // child lines for the - // value of the mother - // line (note that we - // have to multiply by - // two, since the shape - // functions scale with - // the inverse - // Jacobian). we thus - // throw away the - // information of one - // of the child lines, - // but there seems to - // be no other way than - // that... - // - // note: to make things - // consistent, and - // restriction - // independent of the - // order in which we - // travel across the - // cells of the coarse - // grid, we have to - // make sure that we - // take the same small - // line when visiting - // its two neighbors, - // to get the value for - // the mother line. we - // take the first line - // 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.; - this->restriction[0](3,3) = 2.; - - break; - }; - - default: - { - // in case we don't - // have the matrices - // (yet), leave them - // empty. this does not - // prevent the use of - // this FE, but will - // prevent the use of - // these matrices - break; - }; - }; - - break; - }; - - - case 3: // 3d - { - switch (degree) - { - case 1: - { - // same principle as in - // 2d, take one child - // 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.; - this->restriction[3](2,2) = 2.; - - this->restriction[4](4,4) = 2.; - this->restriction[4](7,7) = 2.; - this->restriction[5](5,5) = 2.; - this->restriction[7](6,6) = 2.; - - this->restriction[0](8,8) = 2.; - this->restriction[1](9,9) = 2.; - this->restriction[2](10,10) = 2.; - this->restriction[3](11,11) = 2.; - - break; - }; - - default: - { - // in case we don't - // have the matrices - // (yet), leave them - // empty. this does not - // prevent the use of - // this FE, but will - // prevent the use of - // these matrices - break; - }; - }; - - break; - }; - - default: - Assert (false,ExcNotImplemented()); - } + // 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 @@ -708,6 +502,236 @@ FE_Nedelec<3>::shape_grad_grad_component (const unsigned int i, +template +void +FE_Nedelec::initialize_constraints () +{ + // copy constraint matrices if they + // are defined. otherwise leave + // them at zero size + if (degreeinterface_constraints. + TableBase<2,double>::reinit (this->interface_constraints_size()); + this->interface_constraints.fill (Matrices::constraint_matrices[degree-1]); + }; +} + + + +template +void +FE_Nedelec::initialize_embedding () +{ + if ((degree < Matrices::n_embedding_matrices+1) && + (Matrices::embedding[degree-1][0] != 0)) + 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 + // elements, the row sum must + // be 1, but here the shape + // function is multiplied by + // the inverse of the + // Jacobian, which introduces + // a factor of 1/2 when going + // from mother to child) + for (unsigned int row=0; rowdofs_per_cell; ++row) + { + double sum = 0; + for (unsigned int col=0; coldofs_per_cell; ++col) + sum += this->prolongation[c](row,col); + Assert (std::fabs(sum-.5) < 1e-14, + ExcInternalError()); + }; + }; +} + + + +template +void +FE_Nedelec::initialize_restriction () +{ + switch (dim) + { + case 2: // 2d + { + switch (degree) + { + case 1: + { + // this is a strange + // element, since it is + // both additive and + // then it is also + // not. ideally, we + // would like to have + // the value of the + // shape function on + // the coarse line to + // be the mean value of + // that on the two + // child ones. thus, + // one should make it + // additive. however, + // additivity only + // works if an element + // does not have any + // continuity + // requirements, since + // otherwise degrees of + // freedom are shared + // between adjacent + // elements, and when + // we make the element + // additive, that would + // mean that we end up + // adding up + // contributions not + // only from the child + // cells of this cell, + // but also from the + // child cells of the + // neighbor, and since + // we cannot know + // whether there even + // exists a neighbor we + // cannot simply make + // the element + // additive. + // + // so, until someone + // comes along with a + // better alternative, + // we do the following: + // make the element + // non-additive, and + // simply pick the + // value of one of the + // child lines for the + // value of the mother + // line (note that we + // have to multiply by + // two, since the shape + // functions scale with + // the inverse + // Jacobian). we thus + // throw away the + // information of one + // of the child lines, + // but there seems to + // be no other way than + // that... + // + // note: to make things + // consistent, and + // restriction + // independent of the + // order in which we + // travel across the + // cells of the coarse + // grid, we have to + // make sure that we + // take the same small + // line when visiting + // its two neighbors, + // to get the value for + // the mother line. we + // take the first line + // 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.; + this->restriction[0](3,3) = 2.; + + break; + }; + + default: + { + // in case we don't + // have the matrices + // (yet), leave them + // empty. this does not + // prevent the use of + // this FE, but will + // prevent the use of + // these matrices + break; + }; + }; + + break; + }; + + + case 3: // 3d + { + switch (degree) + { + case 1: + { + // same principle as in + // 2d, take one child + // 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.; + this->restriction[3](2,2) = 2.; + + this->restriction[4](4,4) = 2.; + this->restriction[4](7,7) = 2.; + this->restriction[5](5,5) = 2.; + this->restriction[7](6,6) = 2.; + + this->restriction[0](8,8) = 2.; + this->restriction[1](9,9) = 2.; + this->restriction[2](10,10) = 2.; + this->restriction[3](11,11) = 2.; + + break; + }; + + default: + { + // in case we don't + // have the matrices + // (yet), leave them + // empty. this does not + // prevent the use of + // this FE, but will + // prevent the use of + // these matrices + break; + }; + }; + + break; + }; + + default: + Assert (false,ExcNotImplemented()); + } +} + + + template void FE_Nedelec::initialize_unit_support_points () { diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 8a85752502..aa166ffcc4 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -522,7 +522,18 @@ FE_Q::initialize_embedding () 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 1 + for (unsigned int row=0; rowdofs_per_cell; ++row) + { + double sum = 0; + for (unsigned int col=0; coldofs_per_cell; ++col) + sum += this->prolongation[c](row,col); + Assert (std::fabs(sum-1.) < 1e-14, + ExcInternalError()); + }; + }; } diff --git a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc index 30980f2399..a5106aa7c3 100644 --- a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc +++ b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc @@ -24,6 +24,18 @@ #include +namespace +{ + std::vector + invert_numbering (const std::vector &in) + { + std::vector out (in.size()); + for (unsigned int i=0; i @@ -34,284 +46,44 @@ FE_Q_Hierarchical::FE_Q_Hierarchical (const unsigned int degree) false), std::vector >(FiniteElementData(get_dpo_vector(degree),1).dofs_per_cell, std::vector(1,true))), - degree(degree), - renumber(this->dofs_per_cell, 0), - renumber_inverse(this->dofs_per_cell, 0), - face_renumber(this->dofs_per_face, 0), - polynomial_space(Polynomials::Hierarchical::generate_complete_basis(degree)) + degree(degree), + renumber(lexicographic_to_hierarchic_numbering (*this, degree)), + renumber_inverse(invert_numbering(renumber)), + face_renumber(face_lexicographic_to_hierarchic_numbering (degree)), + polynomial_space(Polynomials::Hierarchical::generate_complete_basis(degree)) { - // do some internal book-keeping on - // cells and faces. if in 1d, the - // face function is empty - lexicographic_to_hierarchic_numbering (*this, degree, renumber); - face_lexicographic_to_hierarchic_numbering (degree, face_renumber); - - // build inverse of renumbering - // vector - for (unsigned int i=0; idofs_per_cell; ++i) - renumber_inverse[renumber[i]]=i; - - // build the dofs_subcell, dofs_cell - // matrices for use in building prolongation, - // restriction, and constraint matrices. - const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line; - - for (unsigned int c=0; c::children_per_cell; ++c) - { - dofs_cell.push_back (FullMatrix (dofs_1d,dofs_1d) ); - dofs_subcell.push_back (FullMatrix (dofs_1d,dofs_1d) ); - - for (unsigned int j=0; j=2)) - { - if (((c==0) && (j==1) && ((k % 2)==0)) || - ((c==1) && (j==0) && ((k % 2)==0))) - dofs_subcell[c](j,k) = -1.; - } - // lower diagonal block - else if ((j>=2) && (k>=2) && (j<=k)) - { - double factor = 1.; - for (unsigned int i=1; i<=j;++i) - factor *= ((double) (k-i+1))/((double) i); - if (c==0) - { - dofs_subcell[c](j,k) = ((k+j) % 2 == 0) ? - std::pow(.5,static_cast(k))*factor : - -std::pow(.5,static_cast(k))*factor; - dofs_cell[c](j,k) = std::pow(2.,static_cast(j))*factor; - } - else - { - dofs_subcell[c](j,k) = std::pow(.5,static_cast(k))*factor; - dofs_cell[c](j,k) = ((k+j) % 2 == 0) ? - std::pow(2.,static_cast(j))*factor : - -std::pow(2.,static_cast(j))*factor; - } - } - } - } - } - // fill constraint matrices - if (dim==2 || dim==3) - { - this->interface_constraints.reinit ( (dim==2) ? 1 + 2*(degree-1) : - 5 + 12*(degree-1) + 4*(degree-1)*(degree-1), - (dim==2) ? (degree+1) : - (degree+1)*(degree+1) ); - switch (dim) - { - case 2: - // vertex node - for (unsigned int i=0; iinterface_constraints(0,i) = dofs_subcell[0](1,i); - // edge nodes - for (unsigned int c=0; c::children_per_cell; ++c) - for (unsigned int i=0; iinterface_constraints(1 + c*(degree-1) + j - 2,i) = - dofs_subcell[c](j,i); - break; - case 3: - for (unsigned int i=0; iinterface_constraints(0,face_renumber[i]) = - dofs_subcell[0](1,i % dofs_1d) * - dofs_subcell[0](1,(i - (i % dofs_1d)) / dofs_1d); - - // boundary vertex nodes - this->interface_constraints(1,face_renumber[i]) = - dofs_subcell[0](1, i % dofs_1d) * - dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(2,face_renumber[i]) = - dofs_subcell[1](1, i % dofs_1d) * - dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(3,face_renumber[i]) = - dofs_subcell[1](0, i % dofs_1d) * - dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(4,face_renumber[i]) = - dofs_subcell[0](0, i % dofs_1d) * - dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); - - // interior edges - for (unsigned int j=0; j<(degree-1); j++) - { - this->interface_constraints(5 + j,face_renumber[i]) = - dofs_subcell[0](1, i % dofs_1d) * - dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + (degree-1) + j,face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * - dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + 2*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](1,i % dofs_1d) * - dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + 3*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](2 + j,i % dofs_1d) * - dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d); - } - - // boundary edges - for (unsigned int j=0; j<(degree-1); j++) - { - // bottom edge - this->interface_constraints(5 + 4*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](2 + j, i % dofs_1d) * - dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + 4*(degree-1) + (degree-1) + j,face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * - dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); - // right edge - this->interface_constraints(5 + 4*(degree-1) + 2*(degree-1) + j,face_renumber[i]) = - dofs_subcell[1](1, i % dofs_1d) * - dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + 4*(degree-1) + 3*(degree-1) + j,face_renumber[i]) = - dofs_subcell[1](1, i % dofs_1d) * - dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); - // top edge - this->interface_constraints(5 + 4*(degree-1) + 4*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](2 + j, i % dofs_1d) * - dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + 4*(degree-1) + 5*(degree-1) + j,face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * - dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); - // left edge - this->interface_constraints(5 + 4*(degree-1) + 6*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](0, i % dofs_1d) * - dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + 4*(degree-1) + 7*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](0, i % dofs_1d) * - dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); - } - - // interior faces - for (unsigned int j=0; j<(degree-1); j++) - for (unsigned int k=0; k<(degree-1); k++) - { - // subcell 0 - this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1),face_renumber[i]) = - dofs_subcell[0](2 + j, i % dofs_1d) * - dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d); - // subcell 1 - this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + (degree-1)*(degree-1),face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * - dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d); - // subcell 2 - this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 2*(degree-1)*(degree-1),face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * - dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d); - // subcell 3 - this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 3*(degree-1)*(degree-1),face_renumber[i]) = - dofs_subcell[0](2 + j, i % dofs_1d) * - dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d); - } - } - break; - } - }; - - // fill prolongation and restriction - // matrices - if (dim==1) - { - for (unsigned int c=0; c::children_per_cell; ++c) - { - this->prolongation[c].reinit (this->dofs_per_cell,this->dofs_per_cell); - this->prolongation[c].fill (dofs_subcell[c]); + // The matrix @p{dofs_cell} contains the + // values of the linear functionals of + // the master 1d cell applied to the + // shape functions of the two 1d subcells. + // The matrix @p{dofs_subcell} constains + // the values of the linear functionals + // on each 1d subcell applied to the + // shape functions on the master 1d + // subcell. + // We use @p{dofs_cell} and + // @p{dofs_subcell} to compute the + // @p{prolongation}, @p{restriction} and + // @p{interface_constraints} matrices + // for all dimensions. + std::vector > + dofs_cell (GeometryInfo<1>::children_per_cell, + FullMatrix (2*this->dofs_per_vertex + this->dofs_per_line, + 2*this->dofs_per_vertex + this->dofs_per_line)); + std::vector > + dofs_subcell (GeometryInfo<1>::children_per_cell, + FullMatrix (2*this->dofs_per_vertex + this->dofs_per_line, + 2*this->dofs_per_vertex + this->dofs_per_line)); + // build these fields, as they are + // needed as auxiliary fields later + // on + build_dofs_cell (dofs_cell, dofs_subcell); + + // then use them to initialize + // other fields + initialize_constraints (dofs_subcell); + initialize_embedding_and_restriction (dofs_cell, dofs_subcell); - this->restriction[c].reinit (this->dofs_per_cell,this->dofs_per_cell); - this->restriction[c].fill (dofs_cell[c]); - } - } - else if (dim==2 || dim==3) - { - for (unsigned int c=0; c::children_per_cell; ++c) - { - this->prolongation[c].reinit (this->dofs_per_cell,this->dofs_per_cell); - this->restriction[c].reinit (this->dofs_per_cell,this->dofs_per_cell); - } - // j loops over dofs in the subcell. - // These are the rows in the - // embedding matrix. - for (unsigned int j=0; jdofs_per_cell; ++j) - { - // i loops over the dofs in the master - // cell. These are the columns in - // the embedding matrix. - for (unsigned int i=0; idofs_per_cell; ++i) - { - switch (dim) - { - case 2: - for (unsigned int c=0; c::children_per_cell; ++c) - { - unsigned int c0 = ((c==1) || (c==2)) ? 1 : 0; - unsigned int c1 = ((c==2) || (c==3)) ? 1 : 0; - - this->prolongation[c](j,i) = - dofs_subcell[c0](renumber_inverse[j] % dofs_1d, - renumber_inverse[i] % dofs_1d) * - dofs_subcell[c1]((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d, - (renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d); - - this->restriction[c](j,i) = - dofs_cell[c0](renumber_inverse[j] % dofs_1d, - renumber_inverse[i] % dofs_1d) * - dofs_cell[c1]((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d, - (renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d); - } - break; - - case 3: - for (unsigned int c=0; c::children_per_cell; ++c) - { - unsigned int c0 = ((c==1) || (c==2) || (c==5) || (c==6)) ? 1 : 0; - unsigned int c1 = ((c==4) || (c==5) || (c==6) || (c==7)) ? 1 : 0; - unsigned int c2 = ((c==2) || (c==3) || (c==6) || (c==7)) ? 1 : 0; - - this->prolongation[c](j,i) = - dofs_subcell[c0](renumber_inverse[j] % dofs_1d, - renumber_inverse[i] % dofs_1d) * - dofs_subcell[c1](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d) % dofs_1d, - ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d) % dofs_1d) * - dofs_subcell[c2](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d - (((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d, - ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d - (((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d); - - this->restriction[c](j,i) = - dofs_cell[c0](renumber_inverse[j] % dofs_1d, - renumber_inverse[i] % dofs_1d) * - dofs_cell[c1](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d) % dofs_1d, - ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d) % dofs_1d) * - dofs_cell[c2](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d - (((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d, - ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d - (((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d); - } - break; - } - } - } - } - else - Assert (false, ExcNotImplemented()); // finally fill in support points // on cell and face initialize_unit_support_points (); @@ -405,6 +177,298 @@ FE_Q_Hierarchical::shape_grad_grad_component (const unsigned int i, //---------------------------------------------------------------------- +template +void +FE_Q_Hierarchical::build_dofs_cell (std::vector > &dofs_cell, + std::vector > &dofs_subcell) const +{ + const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line; + + for (unsigned int c=0; c::children_per_cell; ++c) + for (unsigned int j=0; j=2)) + { + if (((c==0) && (j==1) && ((k % 2)==0)) || + ((c==1) && (j==0) && ((k % 2)==0))) + dofs_subcell[c](j,k) = -1.; + } + // lower diagonal block + else if ((j>=2) && (k>=2) && (j<=k)) + { + double factor = 1.; + for (unsigned int i=1; i<=j;++i) + factor *= ((double) (k-i+1))/((double) i); + if (c==0) + { + dofs_subcell[c](j,k) = ((k+j) % 2 == 0) ? + std::pow(.5,static_cast(k))*factor : + -std::pow(.5,static_cast(k))*factor; + dofs_cell[c](j,k) = std::pow(2.,static_cast(j))*factor; + } + else + { + dofs_subcell[c](j,k) = std::pow(.5,static_cast(k))*factor; + dofs_cell[c](j,k) = ((k+j) % 2 == 0) ? + std::pow(2.,static_cast(j))*factor : + -std::pow(2.,static_cast(j))*factor; + } + } + } +} + + + +template +void +FE_Q_Hierarchical:: +initialize_constraints (const std::vector > &dofs_subcell) +{ + const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line; + + this->interface_constraints + .TableBase<2,double>::reinit (this->interface_constraints_size()); + + switch (dim) + { + case 1: + { + // no constraints in 1d + break; + } + + case 2: + { + // vertex node + for (unsigned int i=0; iinterface_constraints(0,i) = dofs_subcell[0](1,i); + // edge nodes + for (unsigned int c=0; c::children_per_cell; ++c) + for (unsigned int i=0; iinterface_constraints(1 + c*(degree-1) + j - 2,i) = + dofs_subcell[c](j,i); + break; + } + + case 3: + { + for (unsigned int i=0; iinterface_constraints(0,face_renumber[i]) = + dofs_subcell[0](1,i % dofs_1d) * + dofs_subcell[0](1,(i - (i % dofs_1d)) / dofs_1d); + + // boundary vertex nodes + this->interface_constraints(1,face_renumber[i]) = + dofs_subcell[0](1, i % dofs_1d) * + dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(2,face_renumber[i]) = + dofs_subcell[1](1, i % dofs_1d) * + dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(3,face_renumber[i]) = + dofs_subcell[1](0, i % dofs_1d) * + dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(4,face_renumber[i]) = + dofs_subcell[0](0, i % dofs_1d) * + dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); + + // interior edges + for (unsigned int j=0; j<(degree-1); j++) + { + this->interface_constraints(5 + j,face_renumber[i]) = + dofs_subcell[0](1, i % dofs_1d) * + dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(5 + (degree-1) + j,face_renumber[i]) = + dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(5 + 2*(degree-1) + j,face_renumber[i]) = + dofs_subcell[0](1,i % dofs_1d) * + dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(5 + 3*(degree-1) + j,face_renumber[i]) = + dofs_subcell[0](2 + j,i % dofs_1d) * + dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d); + } + + // boundary edges + for (unsigned int j=0; j<(degree-1); j++) + { + // bottom edge + this->interface_constraints(5 + 4*(degree-1) + j,face_renumber[i]) = + dofs_subcell[0](2 + j, i % dofs_1d) * + dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(5 + 4*(degree-1) + (degree-1) + j,face_renumber[i]) = + dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); + // right edge + this->interface_constraints(5 + 4*(degree-1) + 2*(degree-1) + j,face_renumber[i]) = + dofs_subcell[1](1, i % dofs_1d) * + dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(5 + 4*(degree-1) + 3*(degree-1) + j,face_renumber[i]) = + dofs_subcell[1](1, i % dofs_1d) * + dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); + // top edge + this->interface_constraints(5 + 4*(degree-1) + 4*(degree-1) + j,face_renumber[i]) = + dofs_subcell[0](2 + j, i % dofs_1d) * + dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(5 + 4*(degree-1) + 5*(degree-1) + j,face_renumber[i]) = + dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); + // left edge + this->interface_constraints(5 + 4*(degree-1) + 6*(degree-1) + j,face_renumber[i]) = + dofs_subcell[0](0, i % dofs_1d) * + dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); + this->interface_constraints(5 + 4*(degree-1) + 7*(degree-1) + j,face_renumber[i]) = + dofs_subcell[0](0, i % dofs_1d) * + dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); + } + + // interior faces + for (unsigned int j=0; j<(degree-1); j++) + for (unsigned int k=0; k<(degree-1); k++) + { + // subcell 0 + this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1),face_renumber[i]) = + dofs_subcell[0](2 + j, i % dofs_1d) * + dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d); + // subcell 1 + this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + (degree-1)*(degree-1),face_renumber[i]) = + dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d); + // subcell 2 + this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 2*(degree-1)*(degree-1),face_renumber[i]) = + dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d); + // subcell 3 + this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 3*(degree-1)*(degree-1),face_renumber[i]) = + dofs_subcell[0](2 + j, i % dofs_1d) * + dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d); + } + } + break; + } + + default: + Assert (false, ExcNotImplemented()); + } +} + + + +template +void +FE_Q_Hierarchical:: +initialize_embedding_and_restriction (const std::vector > &dofs_cell, + const std::vector > &dofs_subcell) +{ + const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line; + + for (unsigned int c=0; c::children_per_cell; ++c) + { + this->prolongation[c].reinit (this->dofs_per_cell, this->dofs_per_cell); + this->restriction[c].reinit (this->dofs_per_cell, this->dofs_per_cell); + } + + // the 1d case is particularly + // simple, so special case it: + if (dim==1) + { + for (unsigned int c=0; c::children_per_cell; ++c) + { + this->prolongation[c].fill (dofs_subcell[c]); + this->restriction[c].fill (dofs_cell[c]); + } + return; + } + + // for higher dimensions, things + // are a little more tricky: + + // j loops over dofs in the + // subcell. These are the rows in + // the embedding matrix. + // + // i loops over the dofs in the + // master cell. These are the + // columns in the embedding matrix. + for (unsigned int j=0; jdofs_per_cell; ++j) + for (unsigned int i=0; idofs_per_cell; ++i) + switch (dim) + { + case 2: + { + for (unsigned int c=0; c::children_per_cell; ++c) + { + unsigned int c0 = ((c==1) || (c==2)) ? 1 : 0; + unsigned int c1 = ((c==2) || (c==3)) ? 1 : 0; + + this->prolongation[c](j,i) = + dofs_subcell[c0](renumber_inverse[j] % dofs_1d, + renumber_inverse[i] % dofs_1d) * + dofs_subcell[c1]((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d, + (renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d); + + this->restriction[c](j,i) = + dofs_cell[c0](renumber_inverse[j] % dofs_1d, + renumber_inverse[i] % dofs_1d) * + dofs_cell[c1]((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d, + (renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d); + } + break; + } + + case 3: + { + for (unsigned int c=0; c::children_per_cell; ++c) + { + unsigned int c0 = ((c==1) || (c==2) || (c==5) || (c==6)) ? 1 : 0; + unsigned int c1 = ((c==4) || (c==5) || (c==6) || (c==7)) ? 1 : 0; + unsigned int c2 = ((c==2) || (c==3) || (c==6) || (c==7)) ? 1 : 0; + + this->prolongation[c](j,i) = + dofs_subcell[c0](renumber_inverse[j] % dofs_1d, + renumber_inverse[i] % dofs_1d) * + dofs_subcell[c1](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d) % dofs_1d, + ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d) % dofs_1d) * + dofs_subcell[c2](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d - (((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d, + ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d - (((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d); + + this->restriction[c](j,i) = + dofs_cell[c0](renumber_inverse[j] % dofs_1d, + renumber_inverse[i] % dofs_1d) * + dofs_cell[c1](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d) % dofs_1d, + ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d) % dofs_1d) * + dofs_cell[c2](((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d - (((renumber_inverse[j] - (renumber_inverse[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d, + ((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d - (((renumber_inverse[i] - (renumber_inverse[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d); + } + break; + } + + default: + Assert (false, ExcNotImplemented()); + } +} + + template void FE_Q_Hierarchical::initialize_unit_support_points () @@ -548,12 +612,13 @@ FE_Q_Hierarchical::get_dpo_vector(const unsigned int deg) template -void -FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering ( - const FiniteElementData &fe_data, - const unsigned int degree, - std::vector &renumber) +std::vector +FE_Q_Hierarchical:: +lexicographic_to_hierarchic_numbering (const FiniteElementData &fe_data, + const unsigned int degree) { + std::vector renumber (fe_data.dofs_per_cell, 0); + const unsigned int n = degree+1; @@ -624,7 +689,8 @@ FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering ( default: Assert(false, ExcNotImplemented()); } - + + Assert (index < renumber.size(), ExcInternalError()); renumber[index] = i; } }; @@ -709,6 +775,7 @@ FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering ( for (unsigned int jx = 2; jx<=degree ;++jx) { unsigned int tensorindex = tensorstart + jx * incr; + Assert (tensorindex < renumber.size(), ExcInternalError()); renumber[tensorindex] = index++; } } @@ -758,6 +825,7 @@ FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering ( { unsigned int tensorindex = tensorstart + jx * incx + jy * incy; + Assert (tensorindex < renumber.size(), ExcInternalError()); renumber[tensorindex] = index++; } } @@ -772,34 +840,36 @@ FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering ( for (unsigned int jx = 2; jx<=degree; jx++) { const unsigned int tensorindex = jx + jy * n + jz * n * n; + Assert (tensorindex < renumber.size(), ExcInternalError()); renumber[tensorindex]=index++; } } } + + return renumber; } template -void -FE_Q_Hierarchical::face_lexicographic_to_hierarchic_numbering ( - const unsigned int degree, - std::vector &numbering) +std::vector +FE_Q_Hierarchical:: +face_lexicographic_to_hierarchic_numbering (const unsigned int degree) { FiniteElementData fe_data(FE_Q_Hierarchical::get_dpo_vector(degree),1); - FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering (fe_data, - degree, - numbering); + return FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering (fe_data, + degree); } #if (deal_II_dimension == 1) template <> -void -FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int, - std::vector&) -{} +std::vector +FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int) +{ + return std::vector (); +} #endif -- 2.39.5