From baaa8281a127c5abf2f8716e1b958dcf20c092db Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 21 Apr 2003 22:13:02 +0000 Subject: [PATCH] Some re-indenting. git-svn-id: https://svn.dealii.org/trunk@7439 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/source/fe/fe_q_hierarchical.cc | 1156 ++++++++--------- 1 file changed, 578 insertions(+), 578 deletions(-) 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 e161c60731..8f1e3bd1b3 100644 --- a/deal.II/deal.II/source/fe/fe_q_hierarchical.cc +++ b/deal.II/deal.II/source/fe/fe_q_hierarchical.cc @@ -34,11 +34,11 @@ 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(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)) { // do some internal book-keeping on // cells and faces. if in 1d, the @@ -56,259 +56,259 @@ FE_Q_Hierarchical::FE_Q_Hierarchical (const unsigned int degree) 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 (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.; - } + if (((c==0) && (j==1)) || ((c==1) && (j==0))) + dofs_subcell[c](j,k) = .5; + else if (((c==0) && (k==0)) || ((c==1) && (k==1))) + dofs_subcell[c](j,k) = 1.; + else + dofs_subcell[c](j,k) = 0.; + } + // upper right block + else if ((j<=1) && (k>=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; - } + 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.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) { - // center vertex node - this->interface_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); + 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); - } + // 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; } - 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]); + 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]); - this->restriction[c].reinit (this->dofs_per_cell,this->dofs_per_cell); - this->restriction[c].fill (dofs_cell[c]); + 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); - } + 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) + for (unsigned int j=0; jdofs_per_cell; ++j) { - case 2: - for (unsigned int c=0; c::children_per_cell; ++c) + // 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) { - 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); + 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; + } } - 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 @@ -379,7 +379,7 @@ FE_Q_Hierarchical::shape_grad_component (const unsigned int i, template Tensor<2,dim> FE_Q_Hierarchical::shape_grad_grad (const unsigned int i, - const Point &p) const + const Point &p) const { Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); return polynomial_space.compute_grad_grad(renumber_inverse[i], p); @@ -390,8 +390,8 @@ FE_Q_Hierarchical::shape_grad_grad (const unsigned int i, template Tensor<2,dim> FE_Q_Hierarchical::shape_grad_grad_component (const unsigned int i, - const Point &p, - const unsigned int component) const + const Point &p, + const unsigned int component) const { Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); Assert (component == 0, ExcIndexRange (component, 0, 1)); @@ -444,33 +444,33 @@ void FE_Q_Hierarchical::initialize_unit_support_points () 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) - { - if (ix==0) - p(0) = 0.; - else if (ix==1) - p(0) = 1.; - else - p(0) = .5; - if (dim>1) - { - if (iy==0) - p(1) = 0.; - else if (iy==1) - p(1) = 1.; - else - p(1) = .5; - } - if (dim>2) { - if (iz==0) - p(2) = 0.; - else if (iz==1) - p(2) = 1.; + if (ix==0) + p(0) = 0.; + else if (ix==1) + p(0) = 1.; else - p(2) = .5; - } - this->unit_support_points[renumber[k++]] = p; - }; + p(0) = .5; + if (dim>1) + { + if (iy==0) + p(1) = 0.; + else if (iy==1) + p(1) = 1.; + else + p(1) = .5; + } + if (dim>2) + { + if (iz==0) + p(2) = 0.; + else if (iz==1) + p(2) = 1.; + else + p(2) = .5; + } + this->unit_support_points[renumber[k++]] = p; + }; } @@ -503,37 +503,37 @@ void FE_Q_Hierarchical::initialize_unit_face_support_points () 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) - { - if (ix==0) - p(0) = 0.; - else if (ix==1) - p(0) = 1.; - else - p(0) = .5; - if (codim>1) - { - if (iy==0) - p(1) = 0.; - else if (iy==1) - p(1) = 1.; - else - p(1) = .5; - } - if (codim>2) { - if (iz==0) - p(2) = 0.; - else if (iz==1) - p(2) = 1.; + if (ix==0) + p(0) = 0.; + else if (ix==1) + p(0) = 1.; else - p(2) = .5; - } - this->unit_face_support_points[face_renumber[k++]] = p; - }; + p(0) = .5; + if (codim>1) + { + if (iy==0) + p(1) = 0.; + else if (iy==1) + p(1) = 1.; + else + p(1) = .5; + } + if (codim>2) + { + if (iz==0) + p(2) = 0.; + else if (iz==1) + p(2) = 1.; + else + p(2) = .5; + } + this->unit_face_support_points[face_renumber[k++]] = p; + }; } - // we use same dpo_vector as FE_Q + // we use same dpo_vector as FE_Q template std::vector FE_Q_Hierarchical::get_dpo_vector(const unsigned int deg) @@ -549,29 +549,29 @@ 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) + const FiniteElementData &fe_data, + const unsigned int degree, + std::vector &renumber) { 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; - }; + { + 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) - { - Assert (fe_data.dofs_per_vertex == 1, ExcInternalError()); - for (unsigned int i=0; i::vertices_per_cell; ++i) { - unsigned int index = 0; + Assert (fe_data.dofs_per_vertex == 1, ExcInternalError()); + for (unsigned int i=0; i::vertices_per_cell; ++i) + { + unsigned int index = 0; // Find indices of vertices. // Unfortunately, somebody // switched the upper corner @@ -591,190 +591,190 @@ FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering ( // with the complaining // person back then when all // began :-) - switch (dim) - { - case 1: - { - const unsigned int values[GeometryInfo<1>::vertices_per_cell] - = { 0, 1 }; - index = values[i]; - break; - }; + switch (dim) + { + case 1: + { + const unsigned int values[GeometryInfo<1>::vertices_per_cell] + = { 0, 1 }; + index = values[i]; + break; + }; - case 2: - { - const unsigned int values[GeometryInfo<2>::vertices_per_cell] - = { 0, 1, n + 1, n }; - index = values[i]; - break; - }; + case 2: + { + const unsigned int values[GeometryInfo<2>::vertices_per_cell] + = { 0, 1, n + 1, n }; + index = values[i]; + break; + }; - case 3: - { - const unsigned int values[GeometryInfo<3>::vertices_per_cell] - = { 0, 1, - n * n + 1, n * n, - n, n + 1, - n * n + n + 1, n * n + n}; - index = values[i]; - break; - }; + case 3: + { + const unsigned int values[GeometryInfo<3>::vertices_per_cell] + = { 0, 1, + n * n + 1, n * n, + n, n + 1, + n * n + n + 1, n * n + n}; + index = values[i]; + break; + }; - default: - Assert(false, ExcNotImplemented()); - } + default: + Assert(false, ExcNotImplemented()); + } - renumber[index] = i; - } - }; + renumber[index] = i; + } + }; // for degree 2 and higher: Lines, // quads, hexes etc also carry // degrees of freedom if (degree > 1) - { - Assert (fe_data.dofs_per_line == degree-1, ExcInternalError()); - Assert ((fe_data.dofs_per_quad == (degree-1)*(degree-1)) || - (dim < 2), ExcInternalError()); - Assert ((fe_data.dofs_per_hex == (degree-1)*(degree-1)*(degree-1)) || - (dim < 3), ExcInternalError()); - - for (unsigned int i=0; i::lines_per_cell; ++i) { - unsigned int index = fe_data.first_line_index - + i*fe_data.dofs_per_line; - unsigned int incr = 0; - unsigned int tensorstart = 0; + Assert (fe_data.dofs_per_line == degree-1, ExcInternalError()); + Assert ((fe_data.dofs_per_quad == (degree-1)*(degree-1)) || + (dim < 2), ExcInternalError()); + Assert ((fe_data.dofs_per_hex == (degree-1)*(degree-1)*(degree-1)) || + (dim < 3), ExcInternalError()); + + for (unsigned int i=0; i::lines_per_cell; ++i) + { + unsigned int index = fe_data.first_line_index + + i*fe_data.dofs_per_line; + unsigned int incr = 0; + unsigned int tensorstart = 0; // This again looks quite // strange because of the odd // numbering scheme. - switch (i+100*dim) - { - // lines in x-direction - case 100: - case 200: case 202: - case 300: case 302: case 304: case 306: - incr = 1; - break; - // lines in y-direction - case 201: case 203: - case 308: case 309: case 310: case 311: - incr = n; - break; - // lines in z-direction - case 301: case 303: case 305: case 307: - incr = n * n; - break; - default: - Assert(false, ExcNotImplemented()); - } - switch (i+100*dim) - { - // x=y=z=0 - case 100: - case 200: case 203: - case 300: case 303: case 308: - tensorstart = 0; - break; - // x=1 y=z=0 - case 201: - case 301: case 309: - tensorstart = 1; - break; - // y=1 x=z=0 - case 202: - case 304: case 307: - tensorstart = n; - break; - // x=z=1 y=0 - case 310: - tensorstart = n * n + 1; - break; - // z=1 x=y=0 - case 302: case 311: - tensorstart = n * n; - break; - // x=y=1 z=0 - case 305: - tensorstart = n + 1; - break; - // y=z=1 x=0 - case 306: - tensorstart = n * n + n; - break; - default: - Assert(false, ExcNotImplemented()); - } + switch (i+100*dim) + { + // lines in x-direction + case 100: + case 200: case 202: + case 300: case 302: case 304: case 306: + incr = 1; + break; + // lines in y-direction + case 201: case 203: + case 308: case 309: case 310: case 311: + incr = n; + break; + // lines in z-direction + case 301: case 303: case 305: case 307: + incr = n * n; + break; + default: + Assert(false, ExcNotImplemented()); + } + switch (i+100*dim) + { + // x=y=z=0 + case 100: + case 200: case 203: + case 300: case 303: case 308: + tensorstart = 0; + break; + // x=1 y=z=0 + case 201: + case 301: case 309: + tensorstart = 1; + break; + // y=1 x=z=0 + case 202: + case 304: case 307: + tensorstart = n; + break; + // x=z=1 y=0 + case 310: + tensorstart = n * n + 1; + break; + // z=1 x=y=0 + case 302: case 311: + tensorstart = n * n; + break; + // x=y=1 z=0 + case 305: + tensorstart = n + 1; + break; + // y=z=1 x=0 + case 306: + tensorstart = n * n + n; + break; + default: + Assert(false, ExcNotImplemented()); + } - for (unsigned int jx = 2; jx<=degree ;++jx) - { - unsigned int tensorindex = tensorstart + jx * incr; - renumber[tensorindex] = index++; - } - } + for (unsigned int jx = 2; jx<=degree ;++jx) + { + unsigned int tensorindex = tensorstart + jx * incr; + renumber[tensorindex] = index++; + } + } - for (int i=0; i(GeometryInfo::quads_per_cell); ++i) - { - unsigned int index = fe_data.first_quad_index+i*fe_data.dofs_per_quad; - unsigned int tensorstart = 0; - unsigned int incx = 0; - unsigned int incy = 0; - switch (i) - { - // z=0 (dim==2), y=0 (dim==3) - case 0: - tensorstart = 0; incx = 1; - if (dim==2) - incy = n; - else - incy = n * n; - break; - // y=1 - case 1: - tensorstart = n; incx = 1; incy = n * n; - break; - // z=0 - case 2: - tensorstart = 0; incx = 1; incy = n; - break; - // x=1 - case 3: - tensorstart = 1; incx = n; incy = n * n; - break; - // z=1 - case 4: - tensorstart = n * n; incx = 1; incy = n; - break; - // x=0 - case 5: - tensorstart = 0; incx = n; incy = n * n; - break; - default: - Assert(false, ExcNotImplemented()); - } - - for (unsigned int jy = 2; jy<=degree; jy++) - for (unsigned int jx = 2; jx<=degree ;++jx) + for (int i=0; i(GeometryInfo::quads_per_cell); ++i) { - unsigned int tensorindex = tensorstart - + jx * incx + jy * incy; - renumber[tensorindex] = index++; + unsigned int index = fe_data.first_quad_index+i*fe_data.dofs_per_quad; + unsigned int tensorstart = 0; + unsigned int incx = 0; + unsigned int incy = 0; + switch (i) + { + // z=0 (dim==2), y=0 (dim==3) + case 0: + tensorstart = 0; incx = 1; + if (dim==2) + incy = n; + else + incy = n * n; + break; + // y=1 + case 1: + tensorstart = n; incx = 1; incy = n * n; + break; + // z=0 + case 2: + tensorstart = 0; incx = 1; incy = n; + break; + // x=1 + case 3: + tensorstart = 1; incx = n; incy = n * n; + break; + // z=1 + case 4: + tensorstart = n * n; incx = 1; incy = n; + break; + // x=0 + case 5: + tensorstart = 0; incx = n; incy = n * n; + break; + default: + Assert(false, ExcNotImplemented()); + } + + for (unsigned int jy = 2; jy<=degree; jy++) + for (unsigned int jx = 2; jx<=degree ;++jx) + { + unsigned int tensorindex = tensorstart + + jx * incx + jy * incy; + renumber[tensorindex] = index++; + } } - } - if (GeometryInfo::hexes_per_cell > 0) - for (int i=0; i(GeometryInfo::hexes_per_cell); ++i) - { - unsigned int index = fe_data.first_hex_index; + if (GeometryInfo::hexes_per_cell > 0) + for (int i=0; i(GeometryInfo::hexes_per_cell); ++i) + { + unsigned int index = fe_data.first_hex_index; - for (unsigned int jz = 2; jz<=degree; jz++) - for (unsigned int jy = 2; jy<=degree; jy++) - for (unsigned int jx = 2; jx<=degree; jx++) - { - const unsigned int tensorindex = jx + jy * n + jz * n * n; - renumber[tensorindex]=index++; - } - } - } + for (unsigned int jz = 2; jz<=degree; jz++) + for (unsigned int jy = 2; jy<=degree; jy++) + for (unsigned int jx = 2; jx<=degree; jx++) + { + const unsigned int tensorindex = jx + jy * n + jz * n * n; + renumber[tensorindex]=index++; + } + } + } } @@ -782,13 +782,13 @@ FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering ( template void FE_Q_Hierarchical::face_lexicographic_to_hierarchic_numbering ( - const unsigned int degree, - std::vector &numbering) + const unsigned int degree, + std::vector &numbering) { FiniteElementData fe_data(FE_Q_Hierarchical::get_dpo_vector(degree),1); FE_Q_Hierarchical::lexicographic_to_hierarchic_numbering (fe_data, - degree, - numbering); + degree, + numbering); } @@ -797,7 +797,7 @@ FE_Q_Hierarchical::face_lexicographic_to_hierarchic_numbering ( template <> void FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int, - std::vector&) + std::vector&) {} #endif @@ -921,12 +921,12 @@ FE_Q_Hierarchical::get_data (const UpdateFlags update_flags, template void FE_Q_Hierarchical::fill_fe_values ( - const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename Mapping::InternalDataBase &fedata, - FEValuesData &data) const + const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const { // convert data object to internal // data for this class. fails with @@ -963,13 +963,13 @@ FE_Q_Hierarchical::fill_fe_values ( template void FE_Q_Hierarchical::fill_fe_face_values ( - const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, - const unsigned int face, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename Mapping::InternalDataBase &fedata, - FEValuesData &data) const + const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const { // convert data object to internal // data for this class. fails with @@ -1011,14 +1011,14 @@ FE_Q_Hierarchical::fill_fe_face_values ( template void FE_Q_Hierarchical::fill_fe_subface_values ( - const Mapping &mapping, - const typename DoFHandler::cell_iterator &cell, - const unsigned int face, - const unsigned int subface, - const Quadrature &quadrature, - typename Mapping::InternalDataBase &mapping_data, - typename Mapping::InternalDataBase &fedata, - FEValuesData &data) const + const Mapping &mapping, + const typename DoFHandler::cell_iterator &cell, + const unsigned int face, + const unsigned int subface, + const Quadrature &quadrature, + typename Mapping::InternalDataBase &mapping_data, + typename Mapping::InternalDataBase &fedata, + FEValuesData &data) const { // convert data object to internal // data for this class. fails with @@ -1138,7 +1138,7 @@ FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, case 2: { static const unsigned int face_vertices[4][2] = - { {0,1},{1,2},{2,3},{0,3} }; + { {0,1},{1,2},{2,3},{0,3} }; return ((face_vertices[face_index][0] == vertex_no) || (face_vertices[face_index][1] == vertex_no)); @@ -1147,8 +1147,8 @@ FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, case 3: { static const unsigned int face_vertices[6][4] = - { {0,1,2,3},{4,5,6,7},{0,1,5,4}, - {1,5,6,2},{3,2,6,7},{0,4,7,3} }; + { {0,1,2,3},{4,5,6,7},{0,1,5,4}, + {1,5,6,2},{3,2,6,7},{0,4,7,3} }; return ((face_vertices[face_index][0] == vertex_no)|| (face_vertices[face_index][1] == vertex_no)|| (face_vertices[face_index][2] == vertex_no)|| @@ -1180,8 +1180,8 @@ FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, // table technique // again static const unsigned int face_lines[6][4] = - { {0,1,2,3},{4,5,6,7},{0,8,9,4}, - {1,9,5,10},{2,10,6,11},{3,8,7,11} }; + { {0,1,2,3},{4,5,6,7},{0,8,9,4}, + {1,9,5,10},{2,10,6,11},{3,8,7,11} }; return ((face_lines[face_index][0] == line_index)|| (face_lines[face_index][1] == line_index)|| (face_lines[face_index][2] == line_index)|| @@ -1240,102 +1240,102 @@ FE_Q_Hierarchical::get_embedding_dofs (const unsigned int sub_degree) const ExcIndexRange(sub_degree, 1, degree)); if (dim==1) - { - std::vector embedding_dofs (sub_degree+1); - for (unsigned int i=0; i<(sub_degree+1); ++i) - embedding_dofs[i] = i; + { + std::vector embedding_dofs (sub_degree+1); + for (unsigned int i=0; i<(sub_degree+1); ++i) + embedding_dofs[i] = i; - return embedding_dofs; - } + return embedding_dofs; + } if (sub_degree==1) - { - std::vector embedding_dofs (GeometryInfo::vertices_per_cell); - for (unsigned int i=0; i::vertices_per_cell; ++i) - embedding_dofs[i] = i; + { + std::vector embedding_dofs (GeometryInfo::vertices_per_cell); + for (unsigned int i=0; i::vertices_per_cell; ++i) + embedding_dofs[i] = i; - return embedding_dofs; - } + return embedding_dofs; + } else if (sub_degree==degree) - { - std::vector embedding_dofs (this->dofs_per_cell); - for (unsigned int i=0; idofs_per_cell; ++i) - embedding_dofs[i] = i; + { + std::vector embedding_dofs (this->dofs_per_cell); + for (unsigned int i=0; idofs_per_cell; ++i) + embedding_dofs[i] = i; - return embedding_dofs; - } + return embedding_dofs; + } if ((dim==2) || (dim==3)) - { - std::vector embedding_dofs ( (dim==2) ? - (sub_degree+1) * (sub_degree+1) : - (sub_degree+1) * (sub_degree+1) * (sub_degree+1) ); - - for (unsigned int i=0; i<( (dim==2) ? - (sub_degree+1) * (sub_degree+1) : - (sub_degree+1) * (sub_degree+1) * (sub_degree+1) ); ++i) { - // vertex - if (i::vertices_per_cell) - embedding_dofs[i] = i; - // line - else if (i<(GeometryInfo::vertices_per_cell + - GeometryInfo::lines_per_cell * (sub_degree-1))) - { - const unsigned int j = (i - GeometryInfo::vertices_per_cell) % - (sub_degree-1); - const unsigned int line = (i - GeometryInfo::vertices_per_cell - j) / (sub_degree-1); + std::vector embedding_dofs ( (dim==2) ? + (sub_degree+1) * (sub_degree+1) : + (sub_degree+1) * (sub_degree+1) * (sub_degree+1) ); - embedding_dofs[i] = GeometryInfo::vertices_per_cell + - line * (degree-1) + j; - } - // quad - else if (i<(GeometryInfo::vertices_per_cell + - GeometryInfo::lines_per_cell * (sub_degree-1)) + - GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1)) - { - const unsigned int j = (i - GeometryInfo::vertices_per_cell - - GeometryInfo::lines_per_cell * (sub_degree-1)) % (sub_degree-1); - const unsigned int k = ( (i - GeometryInfo::vertices_per_cell - - GeometryInfo::lines_per_cell * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1); - const unsigned int face = (i - GeometryInfo::vertices_per_cell - - GeometryInfo::lines_per_cell * (sub_degree-1) - k * (sub_degree-1) - j) / ( (sub_degree-1) * (sub_degree-1) ); - - embedding_dofs[i] = GeometryInfo::vertices_per_cell + - GeometryInfo::lines_per_cell * (degree-1) + - face * (degree-1) * (degree-1) + - k * (degree-1) + j; - } - // hex - else if (i<(GeometryInfo::vertices_per_cell + - GeometryInfo::lines_per_cell * (sub_degree-1)) + - GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) + - GeometryInfo::hexes_per_cell * (sub_degree-1) * (sub_degree-1) * (sub_degree-1)) - { - const unsigned int j = (i - GeometryInfo::vertices_per_cell - - GeometryInfo::lines_per_cell * (sub_degree-1) - - GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) ) % (sub_degree-1); - const unsigned int k = ( (i - GeometryInfo::vertices_per_cell - - GeometryInfo::lines_per_cell * (sub_degree-1) - - GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1); - const unsigned int l = (i - GeometryInfo::vertices_per_cell - - GeometryInfo::lines_per_cell * (sub_degree-1) - - GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j - k * (sub_degree-1)) / ( (sub_degree-1) * (sub_degree-1) ); + for (unsigned int i=0; i<( (dim==2) ? + (sub_degree+1) * (sub_degree+1) : + (sub_degree+1) * (sub_degree+1) * (sub_degree+1) ); ++i) + { + // vertex + if (i::vertices_per_cell) + embedding_dofs[i] = i; + // line + else if (i<(GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (sub_degree-1))) + { + const unsigned int j = (i - GeometryInfo::vertices_per_cell) % + (sub_degree-1); + const unsigned int line = (i - GeometryInfo::vertices_per_cell - j) / (sub_degree-1); + + embedding_dofs[i] = GeometryInfo::vertices_per_cell + + line * (degree-1) + j; + } + // quad + else if (i<(GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (sub_degree-1)) + + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1)) + { + const unsigned int j = (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1)) % (sub_degree-1); + const unsigned int k = ( (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1); + const unsigned int face = (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1) - k * (sub_degree-1) - j) / ( (sub_degree-1) * (sub_degree-1) ); + + embedding_dofs[i] = GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (degree-1) + + face * (degree-1) * (degree-1) + + k * (degree-1) + j; + } + // hex + else if (i<(GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (sub_degree-1)) + + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) + + GeometryInfo::hexes_per_cell * (sub_degree-1) * (sub_degree-1) * (sub_degree-1)) + { + const unsigned int j = (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1) - + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) ) % (sub_degree-1); + const unsigned int k = ( (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1) - + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1); + const unsigned int l = (i - GeometryInfo::vertices_per_cell - + GeometryInfo::lines_per_cell * (sub_degree-1) - + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j - k * (sub_degree-1)) / ( (sub_degree-1) * (sub_degree-1) ); - embedding_dofs[i] = GeometryInfo::vertices_per_cell + - GeometryInfo::lines_per_cell * (degree-1) + - GeometryInfo::quads_per_cell * (degree-1) * (degree-1) + - l * (degree-1) * (degree-1) + k * (degree-1) + j; - } - } + embedding_dofs[i] = GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (degree-1) + + GeometryInfo::quads_per_cell * (degree-1) * (degree-1) + + l * (degree-1) * (degree-1) + k * (degree-1) + j; + } + } - return embedding_dofs; - } + return embedding_dofs; + } else - { - Assert(false, ExcNotImplemented ()); - return std::vector (); - } + { + Assert(false, ExcNotImplemented ()); + return std::vector (); + } } -- 2.39.5