std::vector<bool> (dim,false),
std::vector<std::vector<bool> >(FiniteElementData<dim>(get_dpo_vector(degree),dim).dofs_per_cell,
std::vector<bool>(dim,true))),
- degree(degree)
+ degree(degree)
{
Assert (dim >= 2, ExcNotUsefulInThisDimension());
switch (dim)
{
case 2: // 2d
- {
- switch (degree)
- {
- case 1:
- {
- // DoF on bottom line
- // of coarse cell will
- // be mean value of
- // bottom DoFs on the
- // two adjacent child
- // cells
- this->restriction[0](0,0) = 0.5;
- this->restriction[1](0,0) = 0.5;
- // same for other DoFs
- this->restriction[1](1,1) = 0.5;
- this->restriction[2](1,1) = 0.5;
-
- this->restriction[2](2,2) = 0.5;
- this->restriction[3](2,2) = 0.5;
-
- this->restriction[3](3,3) = 0.5;
- this->restriction[0](3,3) = 0.5;
-
- break;
- };
+ {
+ switch (degree)
+ {
+ case 1:
+ {
+ // DoF on bottom line
+ // of coarse cell will
+ // be mean value of
+ // bottom DoFs on the
+ // two adjacent child
+ // cells
+ this->restriction[0](0,0) = 0.5;
+ this->restriction[1](0,0) = 0.5;
+ // same for other DoFs
+ this->restriction[1](1,1) = 0.5;
+ this->restriction[2](1,1) = 0.5;
+
+ this->restriction[2](2,2) = 0.5;
+ this->restriction[3](2,2) = 0.5;
+
+ this->restriction[3](3,3) = 0.5;
+ this->restriction[0](3,3) = 0.5;
+
+ break;
+ };
- default:
- {
- // in case we don't
- // have the matrices
- // (yet), set them to
- // impossible
- // values. this does
- // not prevent the use
- // of this FE, but will
- // prevent the use of
- // these matrices
- for (unsigned int i=0;
- i<GeometryInfo<dim>::children_per_cell;
- ++i)
- this->restriction[i].reinit(0,0);
- };
- };
+ default:
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ this->restriction[i].reinit(0,0);
+ };
+ };
- break;
- };
+ break;
+ };
case 3: // 3d
- {
- switch (degree)
- {
- case 1:
- {
- // same principle as in
- // 2d
- this->restriction[0](0,0) = 0.5;
- this->restriction[1](0,0) = 0.5;
+ {
+ switch (degree)
+ {
+ case 1:
+ {
+ // same principle as in
+ // 2d
+ this->restriction[0](0,0) = 0.5;
+ this->restriction[1](0,0) = 0.5;
- this->restriction[1](1,1) = 0.5;
- this->restriction[2](1,1) = 0.5;
+ this->restriction[1](1,1) = 0.5;
+ this->restriction[2](1,1) = 0.5;
- this->restriction[2](2,2) = 0.5;
- this->restriction[3](2,2) = 0.5;
+ this->restriction[2](2,2) = 0.5;
+ this->restriction[3](2,2) = 0.5;
- this->restriction[3](3,3) = 0.5;
- this->restriction[0](3,3) = 0.5;
+ this->restriction[3](3,3) = 0.5;
+ this->restriction[0](3,3) = 0.5;
- this->restriction[4](4,4) = 0.5;
- this->restriction[5](4,4) = 0.5;
+ this->restriction[4](4,4) = 0.5;
+ this->restriction[5](4,4) = 0.5;
- this->restriction[5](5,5) = 0.5;
- this->restriction[6](5,5) = 0.5;
+ this->restriction[5](5,5) = 0.5;
+ this->restriction[6](5,5) = 0.5;
- this->restriction[6](6,6) = 0.5;
- this->restriction[7](6,6) = 0.5;
+ this->restriction[6](6,6) = 0.5;
+ this->restriction[7](6,6) = 0.5;
- this->restriction[7](7,7) = 0.5;
- this->restriction[4](7,7) = 0.5;
+ this->restriction[7](7,7) = 0.5;
+ this->restriction[4](7,7) = 0.5;
- this->restriction[1](8,8) = 0.5;
- this->restriction[5](8,8) = 0.5;
+ this->restriction[1](8,8) = 0.5;
+ this->restriction[5](8,8) = 0.5;
- this->restriction[2](9,9) = 0.5;
- this->restriction[6](9,9) = 0.5;
+ this->restriction[2](9,9) = 0.5;
+ this->restriction[6](9,9) = 0.5;
- this->restriction[3](10,10) = 0.5;
- this->restriction[7](10,10) = 0.5;
+ this->restriction[3](10,10) = 0.5;
+ this->restriction[7](10,10) = 0.5;
- this->restriction[0](11,11) = 0.5;
- this->restriction[5](11,11) = 0.5;
+ this->restriction[0](11,11) = 0.5;
+ this->restriction[5](11,11) = 0.5;
- break;
- };
+ break;
+ };
- default:
- {
- // in case we don't
- // have the matrices
- // (yet), set them to
- // impossible
- // values. this does
- // not prevent the use
- // of this FE, but will
- // prevent the use of
- // these matrices
- for (unsigned int i=0;
- i<GeometryInfo<dim>::children_per_cell;
- ++i)
- this->restriction[i].reinit(0,0);
- };
- };
+ default:
+ {
+ // in case we don't
+ // have the matrices
+ // (yet), set them to
+ // impossible
+ // values. this does
+ // not prevent the use
+ // of this FE, but will
+ // prevent the use of
+ // these matrices
+ for (unsigned int i=0;
+ i<GeometryInfo<dim>::children_per_cell;
+ ++i)
+ this->restriction[i].reinit(0,0);
+ };
+ };
- break;
- };
+ break;
+ };
default:
- Assert (false,ExcNotImplemented());
+ Assert (false,ExcNotImplemented());
}
// finally fill in support points
switch (dim)
{
case 2: // 2D
- {
- switch (degree)
- {
- // first order Nedelec
- // elements
- case 1:
- {
- switch (i)
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements
+ case 1:
{
- // (1-y, 0)
- case 0: return (component == 0 ? 1-p(1) : 0);
- // (0,x)
- case 1: return (component == 0 ? 0 : p(0));
- // (y, 0)
- case 2: return (component == 0 ? p(1) : 0);
- // (0, 1-x)
- case 3: return (component == 0 ? 0 : 1-p(0));
-
- // there are
- // only four
- // shape
- // functions!?
- default:
- Assert (false, ExcInternalError());
- return 0;
+ switch (i)
+ {
+ // (1-y, 0)
+ case 0: return (component == 0 ? 1-p(1) : 0);
+ // (0,x)
+ case 1: return (component == 0 ? 0 : p(0));
+ // (y, 0)
+ case 2: return (component == 0 ? p(1) : 0);
+ // (0, 1-x)
+ case 3: return (component == 0 ? 0 : 1-p(0));
+
+ // there are
+ // only four
+ // shape
+ // functions!?
+ default:
+ Assert (false, ExcInternalError());
+ return 0;
+ };
};
- };
- // no other degrees
- // implemented
- default:
- Assert (false, ExcNotImplemented());
- };
- };
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
case 3: // 3D
- {
- switch (degree)
- {
- // first order Nedelec
- // elements
- case 1:
- {
- // note that the
- // degrees of freedom
- // on opposite faces
- // have a common vector
- // direction, so simply
- // that a little. these
- // directions are:
- //
- // for lines 0, 2, 4, 6:
- // (1,0,0)
- // for lines 1, 3, 5, 7:
- // (0,0,1)
- // for lines 8, 9, 10, 11:
- // (0,1,0)
- //
- // thus, sort out all
- // those cases where
- // the component is
- // zero anyway, and
- // only otherwise
- // compute the
- // spatially dependent
- // part which is then
- // also the return
- // value
- if (((i<8) && (((i%2==0) && (component!=0)) ||
- ((i%2==1) && (component!=2)))) ||
- ((i>=8) && (component != 1)))
- return 0;
-
- // now we know that the
- // only non-zero
- // component is
- // requested:
-//[Anna]:checked: OK
- const double x = p(0),
- y = p(1),
- z = p(2);
- switch (i)
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements
+ case 1:
{
- case 0: return (1-y)*(1-z);
- case 2: return (1-y)*z;
- case 1: return x*(1-y);
- case 3: return (1-x)*(1-y);
-
- case 4: return y*(1-z);
- case 6: return y*z;
- case 5: return x*y;
- case 7: return (1-x)*y;
-
- case 8: return (1-x)*(1-z);
- case 9: return x*(1-z);
- case 10: return x*z;
- case 11: return (1-x)*z;
- default:
- Assert (false, ExcInternalError());
+ // note that the
+ // degrees of freedom
+ // on opposite faces
+ // have a common vector
+ // direction, so simply
+ // that a little. these
+ // directions are:
+ //
+ // for lines 0, 2, 4, 6:
+ // (1,0,0)
+ // for lines 1, 3, 5, 7:
+ // (0,0,1)
+ // for lines 8, 9, 10, 11:
+ // (0,1,0)
+ //
+ // thus, sort out all
+ // those cases where
+ // the component is
+ // zero anyway, and
+ // only otherwise
+ // compute the
+ // spatially dependent
+ // part which is then
+ // also the return
+ // value
+ if (((i<8) && (((i%2==0) && (component!=0)) ||
+ ((i%2==1) && (component!=2)))) ||
+ ((i>=8) && (component != 1)))
return 0;
+
+ // now we know that the
+ // only non-zero
+ // component is
+ // requested:
+//[Anna]:checked: OK
+ const double x = p(0),
+ y = p(1),
+ z = p(2);
+ switch (i)
+ {
+ case 0: return (1-y)*(1-z);
+ case 2: return (1-y)*z;
+ case 1: return x*(1-y);
+ case 3: return (1-x)*(1-y);
+
+ case 4: return y*(1-z);
+ case 6: return y*z;
+ case 5: return x*y;
+ case 7: return (1-x)*y;
+
+ case 8: return (1-x)*(1-z);
+ case 9: return x*(1-z);
+ case 10: return x*z;
+ case 11: return (1-x)*z;
+ default:
+ Assert (false, ExcInternalError());
+ return 0;
+ };
};
- };
- // no other degrees
- // implemented
- default:
- Assert (false, ExcNotImplemented());
- };
- };
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
- // presently no other space
- // dimension implemented
+ // presently no other space
+ // dimension implemented
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
};
return 0;
switch (dim)
{
case 2: // 2D
- {
- switch (degree)
- {
- // first order Nedelec
- // elements
- case 1:
- {
- // on the unit cell,
- // the gradients of
- // these shape
- // functions are
- // constant, so we pack
- // them into a table
- // for simpler lookup
- //
- // the format is: first
- // index=shape function
- // number; second
- // index=vector
- // component, thrid
- // index=component
- // within gradient
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements
+ case 1:
+ {
+ // on the unit cell,
+ // the gradients of
+ // these shape
+ // functions are
+ // constant, so we pack
+ // them into a table
+ // for simpler lookup
+ //
+ // the format is: first
+ // index=shape function
+ // number; second
+ // index=vector
+ // component, thrid
+ // index=component
+ // within gradient
//[Anna]: checked:OK
- static const double unit_gradients[4][2][2]
- = { { {0.,-1.}, {0.,0.} },
- { {0.,0.}, {1.,0.} },
- { {0.,+1.}, {0.,0.} },
- { {0.,0.}, {-1.,0.} } };
- return Tensor<1,dim>(unit_gradients[i][component]);
- };
-
- // no other degrees
- // implemented
- default:
- Assert (false, ExcNotImplemented());
- };
- };
+ static const double unit_gradients[4][2][2]
+ = { { {0.,-1.}, {0.,0.} },
+ { {0.,0.}, {1.,0.} },
+ { {0.,+1.}, {0.,0.} },
+ { {0.,0.}, {-1.,0.} } };
+ return Tensor<1,dim>(unit_gradients[i][component]);
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
- // presently no other space
- // dimension implemented
+ // presently no other space
+ // dimension implemented
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
};
return Tensor<1,dim>();
switch (dim)
{
case 2: // 2D
- {
- switch (degree)
- {
- // first order Nedelec
- // elements. their second
- // derivatives on the
- // unit cell are zero
- case 1:
- {
- return Tensor<2,dim>();
- };
-
- // no other degrees
- // implemented
- default:
- Assert (false, ExcNotImplemented());
- };
- };
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements. their second
+ // derivatives on the
+ // unit cell are zero
+ case 1:
+ {
+ return Tensor<2,dim>();
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
case 3: // 3D
- {
- switch (degree)
- {
- // first order Nedelec
- // elements. their second
- // derivatives on the
- // unit cell are zero
- case 1:
- {
- return Tensor<2,dim>();
- };
-
- // no other degrees
- // implemented
- default:
- Assert (false, ExcNotImplemented());
- };
- };
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements. their second
+ // derivatives on the
+ // unit cell are zero
+ case 1:
+ {
+ return Tensor<2,dim>();
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
- // presently no other space
- // dimension implemented
+ // presently no other space
+ // dimension implemented
default:
- Assert (false, ExcNotImplemented());
+ Assert (false, ExcNotImplemented());
};
return Tensor<2,dim>();