false),
std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),1).dofs_per_cell,
std::vector<bool>(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<double>::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<double>::generate_complete_basis(degree))
{
// do some internal book-keeping on
// cells and faces. if in 1d, the
const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line;
for (unsigned int c=0; c<GeometryInfo<1>::children_per_cell; ++c)
- {
- dofs_cell.push_back (FullMatrix<double> (dofs_1d,dofs_1d) );
- dofs_subcell.push_back (FullMatrix<double> (dofs_1d,dofs_1d) );
-
- for (unsigned int j=0; j<dofs_1d; ++j)
{
- for (unsigned int k=0; k<dofs_1d; ++k)
- {
- // upper diagonal block
- if ((j<=1) && (k<=1))
+ dofs_cell.push_back (FullMatrix<double> (dofs_1d,dofs_1d) );
+ dofs_subcell.push_back (FullMatrix<double> (dofs_1d,dofs_1d) );
+
+ for (unsigned int j=0; j<dofs_1d; ++j)
{
- if (((c==0) && (j==0) && (k==0)) ||
- ((c==1) && (j==1) && (k==1)))
- dofs_cell[c](j,k) = 1.;
- else
- dofs_cell[c](j,k) = 0.;
+ for (unsigned int k=0; k<dofs_1d; ++k)
+ {
+ // upper diagonal block
+ if ((j<=1) && (k<=1))
+ {
+ if (((c==0) && (j==0) && (k==0)) ||
+ ((c==1) && (j==1) && (k==1)))
+ dofs_cell[c](j,k) = 1.;
+ else
+ dofs_cell[c](j,k) = 0.;
- 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.;
- }
+ 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<double>(k))*factor :
- -std::pow(.5,static_cast<double>(k))*factor;
- dofs_cell[c](j,k) = std::pow(2.,static_cast<double>(j))*factor;
- }
- else
- {
- dofs_subcell[c](j,k) = std::pow(.5,static_cast<double>(k))*factor;
- dofs_cell[c](j,k) = ((k+j) % 2 == 0) ?
- std::pow(2.,static_cast<double>(j))*factor :
- -std::pow(2.,static_cast<double>(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<double>(k))*factor :
+ -std::pow(.5,static_cast<double>(k))*factor;
+ dofs_cell[c](j,k) = std::pow(2.,static_cast<double>(j))*factor;
+ }
+ else
+ {
+ dofs_subcell[c](j,k) = std::pow(.5,static_cast<double>(k))*factor;
+ dofs_cell[c](j,k) = ((k+j) % 2 == 0) ?
+ std::pow(2.,static_cast<double>(j))*factor :
+ -std::pow(2.,static_cast<double>(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; i<dofs_1d; ++i)
- this->interface_constraints(0,i) = dofs_subcell[0](1,i);
- // edge nodes
- for (unsigned int c=0; c<GeometryInfo<1>::children_per_cell; ++c)
- for (unsigned int i=0; i<dofs_1d; ++i)
- for (unsigned int j=2; j<dofs_1d; ++j)
- this->interface_constraints(1 + c*(degree-1) + j - 2,i) =
- dofs_subcell[c](j,i);
- break;
- case 3:
- for (unsigned int i=0; i<dofs_1d * dofs_1d; i++)
+ 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)
{
- // 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; i<dofs_1d; ++i)
+ this->interface_constraints(0,i) = dofs_subcell[0](1,i);
+ // edge nodes
+ for (unsigned int c=0; c<GeometryInfo<1>::children_per_cell; ++c)
+ for (unsigned int i=0; i<dofs_1d; ++i)
+ for (unsigned int j=2; j<dofs_1d; ++j)
+ this->interface_constraints(1 + c*(degree-1) + j - 2,i) =
+ dofs_subcell[c](j,i);
+ break;
+ case 3:
+ for (unsigned int i=0; i<dofs_1d * dofs_1d; i++)
+ {
+ // 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);
- // 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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<GeometryInfo<dim>::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; j<this->dofs_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; i<this->dofs_per_cell; ++i)
- {
- switch (dim)
+ for (unsigned int j=0; j<this->dofs_per_cell; ++j)
{
- case 2:
- for (unsigned int c=0; c<GeometryInfo<2>::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; i<this->dofs_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<GeometryInfo<2>::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<GeometryInfo<3>::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<GeometryInfo<3>::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
template <int dim>
Tensor<2,dim>
FE_Q_Hierarchical<dim>::shape_grad_grad (const unsigned int i,
- const Point<dim> &p) const
+ const Point<dim> &p) const
{
Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
return polynomial_space.compute_grad_grad(renumber_inverse[i], p);
template <int dim>
Tensor<2,dim>
FE_Q_Hierarchical<dim>::shape_grad_grad_component (const unsigned int i,
- const Point<dim> &p,
- const unsigned int component) const
+ const Point<dim> &p,
+ const unsigned int component) const
{
Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
Assert (component == 0, ExcIndexRange (component, 0, 1));
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;
+ };
}
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 <int dim>
std::vector<unsigned int>
FE_Q_Hierarchical<dim>::get_dpo_vector(const unsigned int deg)
template <int dim>
void
FE_Q_Hierarchical<dim>::lexicographic_to_hierarchic_numbering (
- const FiniteElementData<dim> &fe_data,
- const unsigned int degree,
- std::vector<unsigned int> &renumber)
+ const FiniteElementData<dim> &fe_data,
+ const unsigned int degree,
+ std::vector<unsigned int> &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<GeometryInfo<dim>::vertices_per_cell; ++i)
{
- unsigned int index = 0;
+ Assert (fe_data.dofs_per_vertex == 1, ExcInternalError());
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ {
+ unsigned int index = 0;
// Find indices of vertices.
// Unfortunately, somebody
// switched the upper corner
// 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<GeometryInfo<dim>::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<GeometryInfo<dim>::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<static_cast<signed int>(GeometryInfo<dim>::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<static_cast<signed int>(GeometryInfo<dim>::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<dim>::hexes_per_cell > 0)
- for (int i=0; i<static_cast<signed int>(GeometryInfo<dim>::hexes_per_cell); ++i)
- {
- unsigned int index = fe_data.first_hex_index;
+ if (GeometryInfo<dim>::hexes_per_cell > 0)
+ for (int i=0; i<static_cast<signed int>(GeometryInfo<dim>::hexes_per_cell); ++i)
+ {
+ unsigned int index = fe_data.first_hex_index;
- for (unsigned int jz = 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++;
+ }
+ }
+ }
}
template <int dim>
void
FE_Q_Hierarchical<dim>::face_lexicographic_to_hierarchic_numbering (
- const unsigned int degree,
- std::vector<unsigned int> &numbering)
+ const unsigned int degree,
+ std::vector<unsigned int> &numbering)
{
FiniteElementData<dim-1> fe_data(FE_Q_Hierarchical<dim-1>::get_dpo_vector(degree),1);
FE_Q_Hierarchical<dim-1>::lexicographic_to_hierarchic_numbering (fe_data,
- degree,
- numbering);
+ degree,
+ numbering);
}
template <>
void
FE_Q_Hierarchical<1>::face_lexicographic_to_hierarchic_numbering (const unsigned int,
- std::vector<unsigned int>&)
+ std::vector<unsigned int>&)
{}
#endif
template <int dim>
void
FE_Q_Hierarchical<dim>::fill_fe_values (
- const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const Quadrature<dim> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &data) const
+ const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const Quadrature<dim> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename Mapping<dim>::InternalDataBase &fedata,
+ FEValuesData<dim> &data) const
{
// convert data object to internal
// data for this class. fails with
template <int dim>
void
FE_Q_Hierarchical<dim>::fill_fe_face_values (
- const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face,
- const Quadrature<dim-1> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &data) const
+ const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face,
+ const Quadrature<dim-1> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename Mapping<dim>::InternalDataBase &fedata,
+ FEValuesData<dim> &data) const
{
// convert data object to internal
// data for this class. fails with
template <int dim>
void
FE_Q_Hierarchical<dim>::fill_fe_subface_values (
- const Mapping<dim> &mapping,
- const typename DoFHandler<dim>::cell_iterator &cell,
- const unsigned int face,
- const unsigned int subface,
- const Quadrature<dim-1> &quadrature,
- typename Mapping<dim>::InternalDataBase &mapping_data,
- typename Mapping<dim>::InternalDataBase &fedata,
- FEValuesData<dim> &data) const
+ const Mapping<dim> &mapping,
+ const typename DoFHandler<dim>::cell_iterator &cell,
+ const unsigned int face,
+ const unsigned int subface,
+ const Quadrature<dim-1> &quadrature,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
+ typename Mapping<dim>::InternalDataBase &fedata,
+ FEValuesData<dim> &data) const
{
// convert data object to internal
// data for this class. fails with
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));
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)||
// 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)||
ExcIndexRange(sub_degree, 1, degree));
if (dim==1)
- {
- std::vector<unsigned int> embedding_dofs (sub_degree+1);
- for (unsigned int i=0; i<(sub_degree+1); ++i)
- embedding_dofs[i] = i;
+ {
+ std::vector<unsigned int> 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<unsigned int> embedding_dofs (GeometryInfo<dim>::vertices_per_cell);
- for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
- embedding_dofs[i] = i;
+ {
+ std::vector<unsigned int> embedding_dofs (GeometryInfo<dim>::vertices_per_cell);
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ embedding_dofs[i] = i;
- return embedding_dofs;
- }
+ return embedding_dofs;
+ }
else if (sub_degree==degree)
- {
- std::vector<unsigned int> embedding_dofs (this->dofs_per_cell);
- for (unsigned int i=0; i<this->dofs_per_cell; ++i)
- embedding_dofs[i] = i;
+ {
+ std::vector<unsigned int> embedding_dofs (this->dofs_per_cell);
+ for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+ embedding_dofs[i] = i;
- return embedding_dofs;
- }
+ return embedding_dofs;
+ }
if ((dim==2) || (dim==3))
- {
- std::vector<unsigned int> 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<GeometryInfo<dim>::vertices_per_cell)
- embedding_dofs[i] = i;
- // line
- else if (i<(GeometryInfo<dim>::vertices_per_cell +
- GeometryInfo<dim>::lines_per_cell * (sub_degree-1)))
- {
- const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell) %
- (sub_degree-1);
- const unsigned int line = (i - GeometryInfo<dim>::vertices_per_cell - j) / (sub_degree-1);
+ std::vector<unsigned int> embedding_dofs ( (dim==2) ?
+ (sub_degree+1) * (sub_degree+1) :
+ (sub_degree+1) * (sub_degree+1) * (sub_degree+1) );
- embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell +
- line * (degree-1) + j;
- }
- // quad
- else if (i<(GeometryInfo<dim>::vertices_per_cell +
- GeometryInfo<dim>::lines_per_cell * (sub_degree-1)) +
- GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1))
- {
- const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell -
- GeometryInfo<dim>::lines_per_cell * (sub_degree-1)) % (sub_degree-1);
- const unsigned int k = ( (i - GeometryInfo<dim>::vertices_per_cell -
- GeometryInfo<dim>::lines_per_cell * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1);
- const unsigned int face = (i - GeometryInfo<dim>::vertices_per_cell -
- GeometryInfo<dim>::lines_per_cell * (sub_degree-1) - k * (sub_degree-1) - j) / ( (sub_degree-1) * (sub_degree-1) );
-
- embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell +
- GeometryInfo<dim>::lines_per_cell * (degree-1) +
- face * (degree-1) * (degree-1) +
- k * (degree-1) + j;
- }
- // hex
- else if (i<(GeometryInfo<dim>::vertices_per_cell +
- GeometryInfo<dim>::lines_per_cell * (sub_degree-1)) +
- GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) +
- GeometryInfo<dim>::hexes_per_cell * (sub_degree-1) * (sub_degree-1) * (sub_degree-1))
- {
- const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell -
- GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
- GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) ) % (sub_degree-1);
- const unsigned int k = ( (i - GeometryInfo<dim>::vertices_per_cell -
- GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
- GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1);
- const unsigned int l = (i - GeometryInfo<dim>::vertices_per_cell -
- GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
- GeometryInfo<dim>::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<GeometryInfo<dim>::vertices_per_cell)
+ embedding_dofs[i] = i;
+ // line
+ else if (i<(GeometryInfo<dim>::vertices_per_cell +
+ GeometryInfo<dim>::lines_per_cell * (sub_degree-1)))
+ {
+ const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell) %
+ (sub_degree-1);
+ const unsigned int line = (i - GeometryInfo<dim>::vertices_per_cell - j) / (sub_degree-1);
+
+ embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell +
+ line * (degree-1) + j;
+ }
+ // quad
+ else if (i<(GeometryInfo<dim>::vertices_per_cell +
+ GeometryInfo<dim>::lines_per_cell * (sub_degree-1)) +
+ GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1))
+ {
+ const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell -
+ GeometryInfo<dim>::lines_per_cell * (sub_degree-1)) % (sub_degree-1);
+ const unsigned int k = ( (i - GeometryInfo<dim>::vertices_per_cell -
+ GeometryInfo<dim>::lines_per_cell * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1);
+ const unsigned int face = (i - GeometryInfo<dim>::vertices_per_cell -
+ GeometryInfo<dim>::lines_per_cell * (sub_degree-1) - k * (sub_degree-1) - j) / ( (sub_degree-1) * (sub_degree-1) );
+
+ embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell +
+ GeometryInfo<dim>::lines_per_cell * (degree-1) +
+ face * (degree-1) * (degree-1) +
+ k * (degree-1) + j;
+ }
+ // hex
+ else if (i<(GeometryInfo<dim>::vertices_per_cell +
+ GeometryInfo<dim>::lines_per_cell * (sub_degree-1)) +
+ GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) +
+ GeometryInfo<dim>::hexes_per_cell * (sub_degree-1) * (sub_degree-1) * (sub_degree-1))
+ {
+ const unsigned int j = (i - GeometryInfo<dim>::vertices_per_cell -
+ GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
+ GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) ) % (sub_degree-1);
+ const unsigned int k = ( (i - GeometryInfo<dim>::vertices_per_cell -
+ GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
+ GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1);
+ const unsigned int l = (i - GeometryInfo<dim>::vertices_per_cell -
+ GeometryInfo<dim>::lines_per_cell * (sub_degree-1) -
+ GeometryInfo<dim>::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j - k * (sub_degree-1)) / ( (sub_degree-1) * (sub_degree-1) );
- embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell +
- GeometryInfo<dim>::lines_per_cell * (degree-1) +
- GeometryInfo<dim>::quads_per_cell * (degree-1) * (degree-1) +
- l * (degree-1) * (degree-1) + k * (degree-1) + j;
- }
- }
+ embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell +
+ GeometryInfo<dim>::lines_per_cell * (degree-1) +
+ GeometryInfo<dim>::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<unsigned int> ();
- }
+ {
+ Assert(false, ExcNotImplemented ());
+ return std::vector<unsigned int> ();
+ }
}
projects / dealii-svn.git / commitdiff